Actually, there's no voodoo, "shift" is defined in the function call
to 0L, so no surprise that:
coverage(IRanges(rstart, rend), shift=1L-start, width=end+shift, ...)
gives the same result as
coverage(IRanges(rstart, rend), shift=1L-start, width=end, ...`)
And I think that my correction is not what is expected from this
function, right? Because "my way" you actually lose the
"offset" (these almost 6Mb in my case). So, what you are doing is
actually storing the offset as the last element of the Rle object,
right? Is that on purpose? That would mean that if I want to plot the
coverage properly for my gene (using the genomic coordinate, i.e.
using a bed or wig file), I would have to move this last element as
the first. Wouldn't it be easier to have an "offset" parameter for the
coverage function which being TRUE would put the offset at the right
place (i.e, the first element of the Rle) and FALSE would simply not
add it?
Best,
On 27 May 2009, at 16:40, Nicolas Delhomme wrote:
Hi all,
(Martin, I guess this one is for you :-))
I have a subset of an AlignedRead
class: AlignedRead
length: 1323 reads; width: 36 cycles
chromosome: 2R 2R ... 2R 2R
position: 5866890 5867511 ... 5867007 5867434
strand: - + ... + -
alignQuality: NumericQuality
alignData varLabels: mismatch
When I call the coverage function, I get back a GenomeData obj. Fine.
Formal class 'GenomeData' [package "BSgenome"] with 4 slots
..@ listData :List of 1
.. ..$ 2R:Formal class 'Rle' [package "IRanges"] with 5 slots
.. .. .. ..@ values : int [1:888] 1 2 5 9 12 14 17 20 21
22 ...
.. .. .. ..@ lengths : int [1:888] 7 1 2 1 1 2 1 3 2 2 ...
.. .. .. ..@ elementMetadata: NULL
.. .. .. ..@ elementType : chr "ANYTHING"
.. .. .. ..@ metadata : list()
..@ elementMetadata: NULL
..@ elementType : chr "ANYTHING"
..@ metadata :List of 6
.. ..$ organism : NULL
.. ..$ provider : NULL
.. ..$ providerVersion: NULL
.. ..$ method : chr "coverage,AlignedRead-method"
.. ..$ coords : chr "leftmost"
.. ..$ extend : int 0
What is surprising is when I look at the Rle object I got:
'integer' Rle instance of length 5868503 with 888 runs
Lengths: 7 1 2 1 1 2 1 3 2 2 ...
Values : 1 2 5 9 12 14 17 20 21 22 ...
which tells me that I have a "covered" region of 5868503 bp instead
of the 1675 bp I'm expecting.
This comes from the last set of length/value from the Rle object:
runLength(test$`2R`)[886:888]
[1] 179 36 5866828
runValue(test$`2R`)[886:888]
[1] 0 1 0
I could track it down to the AlignedRead coverage method (l.179 of
methods-AlignedRead.R in svn rev. 39735), where the code is:
coverage(IRanges(rstart, rend), shift=1L-start, width=end-shift, ...)
As shift is negative, the line should be changed to:
coverage(IRanges(rstart, rend), shift=1L-start, width=end+shift, ...)
But this should not work, as "shift" is not defined beforehand. I
would have expected it to crash, but apparently it is simply is
ignored.... (there's some R voodoo going on there...) and simply
results in this call:
coverage(IRanges(rstart, rend), shift=1L-start, width=end, ...)
Unexpected, but anyway, this corrects it:
coverage(IRanges(rstart, rend), shift=1L-start, width=end+1L-
start, ...)
and gives the same results as
coverage(IRanges(rstart, rend), start=start, end=end, ...)
'integer' Rle instance of length 1675 with 887 runs
Lengths: 7 1 2 1 1 2 1 3 2 2 ...
Values : 1 2 5 9 12 14 17 20 21 22 ...
Best,
---------------------------------------------------------------
Nicolas Delhomme
High Throughput Functional Genomics Center
European Molecular Biology Laboratory
Tel: +49 6221 387 8426
Email: [email protected]
Meyerhofstrasse 1 - Postfach 10.2209
69102 Heidelberg, Germany
_______________________________________________
Bioc-sig-sequencing mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing
---------------------------------------------------------------
Nicolas Delhomme
High Throughput Functional Genomics Center
European Molecular Biology Laboratory
Tel: +49 6221 387 8426
Email: [email protected]
Meyerhofstrasse 1 - Postfach 10.2209
69102 Heidelberg, Germany
_______________________________________________
Bioc-sig-sequencing mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing
