bamstats counts unaligned reads as intergenic.

[IanSudbery]

In the bamstats pipeline, intergenic reads are counted using:

bedtools intersect -a %(tagfile)s
    -b %(contextfile)s
    -bed -v | wc -l
    | xargs printf
    >> %(tmpfile)s_1

The problem with this is that if the %(tagfile)s contains unmapped reads, then they will not overlap any context, and will thus be counted as "intergenic". This leads to a situation where "intergenic" can be bigger than "total".

[jscaber]

Hi @IanSudbery , came across this which is now relevant for me, because I am looking at changes in contextstats in response to an experiment.

I have noticed this in the past as well, and would be keen to fix this. One way I guess would be to simply subtract the number of unmapped reads (like times bp), though we don't have that number easily to hand in the mapping pipeline if I am right? (#85)

The other option would be to create an intergenic bed in the genesets pipeline - my best thought for this would be to take the inverse of transcribed area or one or more of the intergenic bed files from our genesets pipeline and get an actual value obtained by a positive experiment - this might not add up to 100% of mapped bases but is better than the current solution.

Let me know if you have any ideas/preferences (you are probably better versed in this than me), otherwise I'll propose something that seems sensible to me.

Best wishes, Jakub