Hi,

I'm about to analyze ATAC-seq data for the first time and my approach atm is:

1. Calling peaks

2. Merge the narrow peaks files into one

3. Convert the merge bed to saf format for input in featureCounts

4. Using featureCounts (Rsubread package) to get counts

5. Use the counts as input to DESeq2

The thing is that I'm stuck at step 4. Since I cant make sense of the output I get. If I try this on only 2 replicate samples, I get the following:

So first I merge the files and convert them to a .saf file. This then looks like this:

I take my 2 BAM files for the 2 replicates and in R:

fc_PE = featureCounts(files=c(BAM1, BAM2), annot.ext=Merged.saf, isPariedEnd=True). 

If I look at the count values, I notice:

Which I thought was weird. Given that my narrowPeak BED is a file with all my called peaks, how can I have peaks without counts in both my samples? So I looked att the featureCount-objects annotations:

Which gives me the chromosomal regions. Then I went into IGV viewer, and this is "Peak_3" och chr3, (i.e. 0 counts): (this is a zoomed-out, with the peak region clearly marked)

Which don't really looks like a 0-count to me?

So, what I'm doing wrong here? There might be something fundamental I've missed in the workflow I'm trying to conduct here? Maybe I'm using featureCounts in the wrong manner?

Any help highly appreciated!

Best,

JB

Hi JB, I could not clearly see the peak names in your IGV figure but it seems your "Peak_3" overlaps with another peak, ie there is another peak in your SAF annotation that overlaps "Peak_3". featureCounts does not count reads that overlap more than one feature by default and this might be the reason why you got zero counts for "Peak_3".

You may set allowMultiOverlap=TRUE to allow featureCounts to assign multi-overlapping reads to all their overlapping features.