Rsubread counting peaks ATACseq
Entering edit mode
Last seen 6.0 years ago


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!



Rsubread atac-seq narrowpeaks featurecounts • 5.0k views
Entering edit mode

Why isPariedEnd=True ?? In fact, dont you want to count the number of Tn5 nicking sites? I would use isPariedEnd=False even if the data going in is paired end ... i.e. treat each read of the pair as a separate Tn5 nick site. Does that make sense?

Entering edit mode
Wei Shi ★ 3.5k
Last seen 11 days ago
Australia/Melbourne/Olivia Newton-John …

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.

Entering edit mode

Hi Wei,

Thank you for your quick reply. You are indeed correct about the cause of my problem. Sorry for the bad resolution in the picture. The two peaks that you see are overlapping and named differently. allowMultiOverlap=TRUE fixes this, but the real underlying problem in my case is of course the merge of the peak files. I want those two very close peaks to be merged as one. Turns out the problem was that I didn't have sorted peak files when using bedops to merge them. 

Your reply made me realize this, thank you!





Login before adding your answer.

Traffic: 486 users visited in the last hour
Help About
Access RSS

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6