Hi,
I am trying to calculate the cut sites for my ATAC-seq data and then use these for featureCounts. However, I am unclear exactly how to write the cut sites in BED file format. I have made two different BED files.
In the first, I modified the start and stop position for all forward (+) reads to be the original start site, and I modified the start and stop of all the reverse (-) reads to be the original stop site - 1. The reason for this is because from my reading, the stop site in a BED file is not inclusive.
BED file 1 e.g.
chr1 10073 10073 NB502048:39:HT3HMBGX3:1:13311:12074:5101/1 30 +
chr1 10073 10073 NB502048:39:HT3HMBGX3:4:21410:17550:14082/2 30 +
chr1 10162 10162 NB502048:39:HT3HMBGX3:2:23103:14543:10456/2 30 +
chr1 10208 10208 NB502048:39:HT3HMBGX3:2:23103:14543:10456/1 30 -
In the second, I modified the start position for all forward reads to be the original start site and the stop position to be start site + 1. For the reverse reads, I modified the start to be the original stop site - 1, and the end to be the original stop site. So in the first BED file, the start and end are the same, whereas in the second, the end is 1 bp more than the start.
BED file 2 e.g.
chr1 10073 10074 NB502048:39:HT3HMBGX3:1:13311:12074:5101/1 30 +
chr1 10073 10074 NB502048:39:HT3HMBGX3:4:21410:17550:14082/2 30 +
chr1 10162 10163 NB502048:39:HT3HMBGX3:2:23103:14543:10456/2 30 +
chr1 10208 10209 NB502048:39:HT3HMBGX3:2:23103:14543:10456/1 30 -
After creating these cut-site BED files, I converted them to BAMs and used these for featureCounts. With the two files, I get slightly different counts. I was expecting them to be the same if the stop site is non-inclusive in BED files.
Could anyone explain which of these approaches (if any) is correct and why I see different results with featureCounts using the two files?
Many thanks for the help,
Lucy
