[Bioc-devel] Confused by readGAlignments(scanBamParam( which ))

[Leonardo Collado Torres]

Hi,

It took me a while, but I was confused by getting different results
when using 'which' in reading BAM files. In my use case, I have a
GRanges with 15 ranges (different exons from the transcripts of a
gene) and I get much higher coverage values than simply using the
range of the 15 ranges.

I now realize that by using all 15 ranges I was duplicating the number
of alignments to use for calculating the coverage instead of what I
wanted: get all the reads that mapped to any of the 15 ranges counting
them once each.

You could say that as a user I was incorrectly using the function. Or
that it would be good to clarify that if you use overlapping ranges in
'which' then you could read a sequence read more than once.

Anyhow, I think that clarifying would be useful to others in the
future. Even a simple example could show this situation so others are
aware.

The code and output is available at
https://gist.github.com/lcolladotor/76f2efdd5d5108772696

Cheers,
Leo