Peak calling software, such as MACS2 or HMMRATAC or Genrich have been developed to identify regions of enrichment and assign statistical metrics to their results. You can choose to believe their assignments or, as you describe, you can choose to not fully believe them and take additional steps to alter the results, such as merging peaks that are close to each other. There is no universal answer to your question (in fact, I would ask first why you would want to merge the peaks, i.e., what are you gaining for your downstream analyses?).

After I posted the original question, I read NucleoATAC paper by Schep AN et al. In nucleoATAC, one of the step includes extending identifying region by "bedtools slop". I extended for 200bp. I get good results now.

Following are links posted previously which I missed.

Nucleoatac - extending broadPeak regions

https://github.com/GreenleafLab/NucleoATAC/issues/46

Thanks for comments.

You can use GenomicRanges in R to merge peaks (as either the union or intersection, depending upon what you are looking for).

For example, let's say you have a .bed file (without a header), and you can create a GenomicRanges object as follows:

gr1 = reduce(GRanges(Rle(input.table$V1),
        IRanges(start=input.table$V2, end=input.table$V3)))

In this situation, I don't think the reduce command is necessary (to merge overlapping rows), but it won't hurt.

If you wanted to then create the union of two peaks, you could use:

peakGr = union(gr1, gr2)

Of, if you wanted to create the intersection of two peaks, you could use:

peakGr = intersect(gr1, gr2)

You can then quantify counts with featureCounts, htseq-count, etc. Of, if you only have 2 samples, you can try comparing one sample as the reference for peak calling (for what is intended to be used for "input" samples).