Hi all,

Recently I am trying to call differential peaks using DEseq2. I have got ATAC-seq data for 10 samples  (5 time points and each time point has 2 replicates). Here is the workflow I trying to establish:

1) Call peaks using MACS2

2) Count the number of reads for each peak in each sample separately using featureCounts (Rsubread package). The output matrix is as follow.

           sample1

peak_1    100 

peak_2   178

3) Merge 10 read counts matrices

           sample1  sample2  sample3……

peak_1    100        94

peak_2   178        243

4) Feed the merge matrix to DEseq2

However, I encounter a trouble when I am trying to merge these 10 matrices Because the rownames of each matrix is different. Some of them represent the same or similar location of the chromosome but using different peak id. I am wondering what's the best way to solve this problem. I will appreciate any hits! Thanks!

The DiffBind package is actively being used for this type of analysis of ATAC-seq data. You can feed in your peaks and reads and it will form the peak matrix as your describe, populate it with read counts, and run a DESeq2 analysis. It also includes useful reporting and plotting functions.

 The accompanying vignette provides clear examples of how to use it.

Cheers-

Rory