Hmm, I actually tried this as a toy example and it works fine. Interval 1 has length 100 because 200-100=100, intersect length is 200-150=50 and subtract length is 150-100=50 which sum up to 100. The issue comes up with larger files, so I'm not sure how to even debug.
Hi,
I have a BED file of features in the genome, for which I am trying to see how many bases overlap exons and how many don't. I used bedtools intersect to get a BED file that gives me the intersection of my features with an exon annotation from gencode, and bedtools sutract to get a BED file for the bases in my features that are not in exons. I then used awk -F'\t' 'BEGIN{SUM=0}{ SUM+=$3-$2 }END{print SUM}' to get the number of bases covered by each of these files but I am running into the strange issue that the bases of the intersection + the length of the subtraction do not add up to the number of bases in the original file. They add up to slightly more than the original.
Am I making some basic logic error, misunderstanding the intersect and subtract functionality or is something very strange going on?
Thanks,
Adriana
1 answer
It's because the endpoint in Bedtools is included in both data sets. For example, if interval 1 is 100-200 and interval 2 is 150-250, bedtools intersect returns 150-200, while subtract returns 100-150. So position 150 gets counted twice when you sum the two data sets.
EDIT: This explanation is incorrect. See OP's reply.
Sorry, you're right. I didn't read your awk command carefully. I was thinking in terms of string lengths.
If your exon annotation files contains isoforms/overlapping sequences, that will cause discrepancies in the numbers. E.g., for a feature interval 100-400 and exon intervals 150-250 and 200-300, intersect returns both 150-250 and 200-300, while subtract returns 100-150 and 300-400. The sum of those intervals is 350, not 300.
Log in to answer this question.
Figured it out. The exons sometimes overlap resulting in the intersection being reported twice. Had to merge the exon annotation and everything worked.