I'm trying to count paired end ATAC-seq alignments that fall within peaks specified by a BED file for differential accessibility analyses across two conditions. Peaks were called using MACS2 with --shift -75 --extsize 150 --nomodel --keep-dup all --SPMR, to center reads around sites of transposase insertion. My question is about applying featureCounts to identify the number of transposition events in each BED-specified interval: I'm not supplying the -p arg to featureCounts because I don't care where the fragments are, but rather where the transposition sites fall. Assuming this is a correct way to think about this, subread v1.6.2 provides a --read2pos 5 option for featureCounts. Should I be using this to specifically count where the transposition events are occurring, rather than the overlap of library fragments over intervals? Alternatively, to keep things consistent with the MACS2 peak calling, is it better to call featureCounts with the same --shiftsize and --extsize parameters I used with MACS2?
