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!