Hi! 
Now I have 6 ATAC-seq samples which have been aligned to hg19 and called peaks by MACS2. 3 of them are from A cells and the other are from B cells. Now I'm going to do differential peak calling by MACS2. Before it, I need to merge peaks from 3 samples of the same cell type. 
I use this command to call peaks.

    macs2 callpeak -B -t sample_1.bam -n sample_1 --nomodel --extsize 120 --shift 60  \
    --outdir sample_1 --keep-dup all

And then I get 5 output files.

    sample_1_control_lambda.bdg
    sample_1_peaks.narrowPeak
    sample_1_peaks.xls
    sample_1_summits.bed
    sample_1_treat_pileup.bdg

But I'm confused how to merge them. 
Shall I use *bedtools intersect* to merge 3 control_lambda.bdg files and  treat_pileup.bdg files of A cell respectively? I'm worried if it will influence the value in the fourth column of the bdg files. 

Or shall I merge their bam files before peak calling?








In my opinion "differential peak calling" does not exist with `macs2` as it is not able to make proper use of experimental replication and as such is not suitable to derive differential calls. It is great at calling peaks because it has been developed for that, but not more than that. I would call peaks on the merged files:

    macs2 callpeak -t *.bam (...)

There is no need to merge the bam files as `t` can read multiple bams and then combine the pipeup internally.

Then convert the narrowPeak file to SAF format:

https://www.biostars.org/p/228636/

Then make a count matrix with `featureCounts` from the subread package:

    featureCounts -a your.saf -F SAF -o out.counts *.bam

...and then feed this into something like DESeq2 or edgeR. Alternatively, look into the [DiffBind][1] or [csaw][2] packages which will instruct you on how to run a proper differential analysis. 


  [1]: https://bioconductor.org/packages/release/bioc/vignettes/DiffBind/inst/doc/DiffBind.pdf
  [2]: http://bioconductor.org/books/3.13/csawBook/