On Wed, Aug 25, 2010 at 12:16 AM, Patrick Aboyoun <paboyoun@fhcrc.org>wrote:

>  Steve,
> Riffing off of Martin's post, I would use the match function to do the
> mapping:
>
> fence.posts <- reduce(reads)
> ranges(reads) <- ranges(fence.posts)[match(reads, fence.posts)]
>
>
This is the first one that came to mind and is the cleanest, but a cleaner
version of the cut method using base R is findInterval:

ranges(reads) <- ranges(fence.posts)[findInterval(start(reads),
start(fence.posts))]

Will be a bit faster than match (which is just findOverlaps(multiple=FALSE))


It is not as fast as the Rle/cut based method by Martin, but it is more
> readable and, thus, easier to maintain.
>
>
> Cheers,
> Patrick
>
>
>
> On 8/24/10 9:17 PM, mtmorgan@fhcrc.org wrote:
>
>> Quoting Steve Lianoglou <mailinglist.honeypot@gmail.com>:
>>
>>  In case anyone else is trying to sharpen their IRanges-fu, here's a
>>> possible solution to my own question -- works much faster.
>>>
>>> Instead of:
>>>
>>>   reads <- GRanges(seqnames='chr1', strand='+',
>>>>               ranges=IRanges(start=c(sample(18:22, 10, replace=TRUE),
>>>>                                sample(100:108, 10, replace=TRUE)),
>>>>                 width=sample(18:20, 20, replace=TRUE)))
>>>>
>>>> You can see there are two "islands" of reads there, at around
>>>> positions 20-40, and 100-120.
>>>> What I want to do is to "fix" the reads in each island to have the
>>>> same start/end position.
>>>>
>>>> Here's code I have to do that (forget any "strand" issues here --
>>>> those are already taken care of):
>>>>
>>>>  fence.posts <- reduce(reads)
>>>>  o <- findOverlaps(fence.posts, reads)
>>>>
>>>
>> Just
>>
>>  ranges(reads)[subjectHits(o)] <- ranges(fence.posts)[queryHits(o)]
>>
>> ? (I would have done findOverlaps(reads, fence.posts), but...)
>>
>>   subjects <- subjectHits(o)
>>>>  query <- queryHits(o)
>>>>
>>>>  for (idx in unique(query)) {
>>>>   adjust <- subjects[query == idx]
>>>>   start(reads[adjust]) <- start(fence.posts[idx])
>>>>   end(reads[adjust]) <- end(fence.posts[idx])
>>>>  }
>>>>
>>>
>>> I'm doing:
>>>
>>> covr <- coverage(reads)[[1]]
>>> fences <- slice(covr, 1)
>>>
>>
>> or
>>
>>  s = covr > 1
>>  breaks = cumsum(runLength(s))[!runValue(s)]
>>  post = cut(start(reads), c(breaks, Inf), labels=FALSE)
>>  ranges(reads) = ranges(fence.posts)[post]
>>
>> ? but the method above keeps track of sequences...
>>
>> Martin
>>
>>  f.repeat <- viewMaxs(fences)
>>> repacked <- IRanges(start=rep(start(fences), f.repeat),
>>>                    end=rep(end(fences), f.repeat))
>>> if (length(repacked) != length(reads)) {
>>>  stop("Length of repacked reads is not the same as the original")
>>> reads <- GRanges(seqnames=seqnames(reads), ranges=repacked,
>>>                 strand=strand(reads))
>>>
>>> I'm assuming that all the reads come from the same chromosome (so I
>>> unlist the result from coverage), and I'm only sending in reads of the
>>> same strand to this part of the function.
>>>
>>> Anyway, running this newer version on 1000 tags takes 0.181 seconds,
>>> where the older version took 20 seconds.
>>>
>>> Thought some might find it interesting, otherwise sorry to spam.
>>>
>>> -steve
>>>
>>> --
>>> Steve Lianoglou
>>> Graduate Student: Computational Systems Biology
>>>  | Memorial Sloan-Kettering Cancer Center
>>>  | Weill Medical College of Cornell University
>>> Contact Info: http://cbio.mskcc.org/~lianos/contact<http://cbio.mskcc.org/%7Elianos/contact>
>>>
>>> _______________________________________________
>>> Bioconductor mailing list
>>> Bioconductor@stat.math.ethz.ch
>>> https://stat.ethz.ch/mailman/listinfo/bioconductor
>>> Search the archives:
>>> http://news.gmane.org/gmane.science.biology.informatics.conductor
>>>
>>>
>> _______________________________________________
>> Bioconductor mailing list
>> Bioconductor@stat.math.ethz.ch
>> https://stat.ethz.ch/mailman/listinfo/bioconductor
>> Search the archives:
>> http://news.gmane.org/gmane.science.biology.informatics.conductor
>>
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor@stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives:
> http://news.gmane.org/gmane.science.biology.informatics.conductor
>

	[[alternative HTML version deleted]]

