Question: Rsubread counting peaks ATACseq
gravatar for j.bergenstrahle
23 months ago by
j.bergenstrahle80 wrote:


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!



ADD COMMENTlink modified 23 months ago by Wei Shi3.2k • written 23 months ago by j.bergenstrahle80

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?

ADD REPLYlink written 20 hours ago by vishal.koparde0
Answer: Rsubread counting peaks ATACseq
gravatar for Wei Shi
23 months ago by
Wei Shi3.2k
Wei Shi3.2k wrote:

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.

ADD COMMENTlink written 23 months ago by Wei Shi3.2k

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!




ADD REPLYlink written 23 months ago by j.bergenstrahle80
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 165 users visited in the last hour