On Mon, Feb 3, 2014 at 4:38 AM, Philip Lijnzaad <p.lijnzaad@umcutrecht.nl>wrote:

> On 2014-02-03 13:07, Philip Lijnzaad wrote:
>
>> Dear all,
>>
>> I can't seem to select a run of scores from GRanges object, using the
>> coordinates from another GRanges object.
>>
>> An example:
>>
>> coverage <- GRanges(seqnames="chrI",
>>
>>              ranges=IRanges(start=c(1,11,21), width=10),
>>
>>              score=c(0, 42, 0)
>>
>>              )
>>
>> feature <- GRanges(seqnames="chrI", ranges=IRanges(start=5, end=25))
>>
>> In ascii this looks like:
>>
>> coverage: ...__________----------__________...
>>
>> feature:  .......xxxxxxxxxxxxxxxxxxxxx........
>>
>> wanted:          _____----------______
>>
>> I tried using
>>
>>    score(subsetByOverlaps(coverage, feature))
>>
>> which I hope  would give me back
>>
>>    0  0  0  0  0 42 42 42 42 42 42 42 42 42 42  0  0  0  0  0
>>
>> or maybe
>>
>>   numeric-Rle of length 20 with 3 runs
>>
>>     Lengths:  5 10  5
>>
>>     Values :  0 42  0
>>
>
> Mmm ... seems I stumbled across a solution myself. If I do
>
>    scores <- Rle(score(coverage), width(ranges(coverage)))
>
>
> I can do
>
>   scores[ ranges(feature) ]
>
>
> which seems to do exactly what I need. Is this the right (and fastest)
> way? If so, might be nice to have this
> in one of the vignettes.
>
>
Note that this solution assumes that the data are on one chromosome. In
general, GRanges maps to RleList, so instead of:

scores <- Rle(score(coverage), width(coverage))
scores[ ranges(feature) ]

You want:

scores <- coverage(coverage, weight="score")
scores[as(feature, "RangesList")]

Make sure that the seqlevels are the same between "coverage" and "feature".

However, in this case fastest way is to have rtracklayer load the coverage
directly as an Rle. In release, this is just

scores <- import("my.bw", asRle=TRUE)
scores[ as(feature, "RangesList") ]

In devel, the first line needs to be

scores <- import("my.bw", as="RleList")

and there is also support for as="NumericList" thanks to Val.

If you are *only* interested in the data over "feature", then you could
tell rtracklayer to only load that data (the rest of the RleList will be
zeroed out).

scores <- import("my.bw", asRle=TRUE, which=feature)
feature.scores <- scores[ as(feature, "RangesList") ]

If "feature" is contained within one chromosome, then you might want to
simply:
unlist(feature.scores, use.names=FALSE)

With all that explained, it would be nice to know your use case. Often
these types of problems are solvable with RleViews[List].

Michael

Cheers,
>
> Ph
>
>
> ------------------------------------------------------------
> ------------------
>
> De informatie opgenomen in dit bericht kan vertrouwelijk zijn en is
> uitsluitend bestemd voor de geadresseerde. Indien u dit bericht onterecht
> ontvangt, wordt u verzocht de inhoud niet te gebruiken en de afzender
> direct
> te informeren door het bericht te retourneren. Het Universitair Medisch
> Centrum Utrecht is een publiekrechtelijke rechtspersoon in de zin van de
> W.H.W.
> (Wet Hoger Onderwijs en Wetenschappelijk Onderzoek) en staat geregistreerd
> bij
> de Kamer van Koophandel voor Midden-Nederland onder nr. 30244197.
>
> Denk s.v.p aan het milieu voor u deze e-mail afdrukt.
>
> ------------------------------------------------------------
> ------------------
>
> This message may contain confidential information and ...{{dropped:13}}

