I have a bed file with regions of interest and a bam file from STAR - aligned RNA-Seq sample. I want to create a small bash script that would result in a txt file that reports position-wise coverage in these regions.
I have been playing with bedtools: I intersected the bed and bam and then I put the result into
bedtools genomecov. The problem is that this tool creates a coverage for whole genome, so even though the bam has been intersected with the regions in bed and I end up with a file that has position-wise coverage it is 95% zeros, obviously...
I think if I use
bedtools coverage this might work but then I noticed a problem described here:
In order to be consistent with IGV I would like to exclude the gaps from alignments.
Any ideas how to tackle this in another way?