Re: [Bioc-sig-seq] A bug in the AlignedRead coverage method?

Martin Morgan Thu, 28 May 2009 09:03:55 -0700

Thanks Nicolas; in some ways I'm tempted to change the error and requirethe user to provide an integer value (append 'L', for instance, oras.integer) to avoid ambiguities or mask user errors through internalconversion. I'll give it some thought.


Martin


Nicolas Delhomme wrote:

Hi Martin,

Thanks, it works nicely. However, when I added a width parameter, itkept on telling me


Error: UserArgumentMismatch
  'names(width)' (or 'names(end)') mismatch with 'levels(chromosome(x))'
  see ?"AlignedRead-class"

whereas my vector has the proper names.

I dig up and the problem comes from the fact that my vector contains anumeric value and not an integer.

In the coverage function, you convert numeric width values into integerones and this removes the names attrs from the vector. (see below for anexample)


 > widths=c(21146708)
 > names(widths)<-"2R"
 > str(widths)
 Named num 21146708
 - attr(*, "names")= chr "2R"
 > as.integer(widths)
[1] 21146708

Adding these two lines fixes it (around l.191 0f methods-AlignedRead.R):

    # save the names (added line)
        w.names<-names(width)

    # convert to integer (this was l.191)
            width <- as.integer(width)

    # restore the names (added line)
        names(width)<-w.names

Hope this helps,

---------------------------------------------------------------
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
---------------------------------------------------------------



On 28 May 2009, at 05:47, Martin Morgan wrote:

Hi Nicolas --

Thanks for your report; S4 is I think evaluating a function argument
in the wrong frame.

In the RELEASE version of ShortRead, v. 1.2.1 available in svn now and
via biocLite in the next 48 hours, the idea is that

 coverage(aln)

trims leading and trailing nucleotides with 0 coverage, so the
coordinates of the Rle are relative to the first covered nucleotide
(in v. 1.2.0 and devel <1.3.8 the tail of the Rle included a run of
unintended 0 coverage). Genomic coordinates, suitable for bed files
etc., are obtained with

 coverage(aln, start=1L)

In the DEVEL version of ShortRead, v. 1.3.9, the default behavior and
preferred way of specifying coordinates has been changed to be
consistent with recent changes in IRanges and Biostrings.

 coverage(aln)

returns genomic (1-based) coordinates, and the prefered way to adjust
the coordinate system is with shift/width arguments (which when
present are named integer vectors, with all levels(chromosome(aln))
present). See ?"AlignedRead-class" and its example for additional
detail.

Thanks again for the report. Please let me know if I've
mis-understood.

Martin

Nicolas Delhomme <[email protected]> writes:

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


--
Martin Morgan
Computational Biology / Fred Hutchinson Cancer Research Center
1100 Fairview Ave. N.
PO Box 19024 Seattle, WA 98109

Location: Arnold Building M1 B861
Phone: (206) 667-2793



--
Martin Morgan
Computational Biology / Fred Hutchinson Cancer Research Center
1100 Fairview Ave. N.
PO Box 19024 Seattle, WA 98109

Location: Arnold Building M1 B861
Phone: (206) 667-2793

_______________________________________________
Bioc-sig-sequencing mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing

Re: [Bioc-sig-seq] A bug in the AlignedRead coverage method?

Reply via email to