Hi Yannick,
On 10/02/2012 10:07 AM, Yannick Wurm wrote:
> Hello,
>
> my Granges object "bigGranges" contains too much data. THat makes
plotting extremely slow. Sometimes I want an averaged overview. I use
the following code below to obtain it. Subsequently I can use ggbio
for plotting.
> But I feel there must be an elegant one-line manner of doing this?
The most elegant manner would be that ggbio takes care of this for
you ;-)
The code you provide below is not self-contained (we don't know what
bigGranges and myscaffolds are). Also you don't explain what it's
supposed to do and/or what the output is supposed to look like.
Anyway, using 2 big nested "for" loops and growing the averageGranges
object at each iteration is in general very inefficient. The "R way"
is
to take advantage of the fact that most arithmetic operations in R are
vectorized. Also the IRanges/GenomicRanges packages provide a lot of
arithmetic and other operations to let you work on RleViews,
RleViewsList, IRangesList, GRanges etc... objects in a vectorized
fashion.
But most importantly, the kind of average you compute seems
questionable
to me. Right now, if a bin is for example the [500, 599] interval, and
if your bigGranges object has 2 ranges that fall within that bin e.g.
start end myValue
510 510 2.8
500 549 3.8
then you will assign the non-weighted myValue average to the bin (i.e.
3.3). Is it really what you want? Or, given that the 1st range doesn't
count a lot, and the 2nd range covers half of the bin, wouldn't you
want the average to be close to 1.9?
Another (somewhat related) issue I see is that if a range in
bigGranges
overlaps with 2 bins then you're counting it twice (because you're
using
subsetByOverlaps()).
What about computing the weighted coverage of bigGranges (using the
myValue metadata column for the weights), and then the mean of this
coverage for a set of bins you specifiy? You would not only get rid
of the above issues, but it can also be implemented in a very
efficient way.
Here is what it would look like:
### First define the bins. The bins would typically be stored in a
### GRanges or IRangesList object. Here we choose the latter:
mybins <- lapply(seqlengths(bigGranges),
function(chrom_len)
IRanges(breakInChunks(chrom_len, 200L)))
mybins <- IRangesList(mybins)
### Then compute the coverage of your 'bigGranges', weighted by
### the values in its "myValue" metadata column:
cvg <- coverage(bigGranges, weight="myValue")
### Then put views (corresponding to the bins) on top of the
### RleList object returned by coverage(). The result is stored
### in a RleViewsList object:
viewslist <- RleViewsList(
lapply(names(cvg),
function(chrom)
Views(cvg[[chrom]], mybins[[chrom]])))
names(viewslist) <- names(cvg)
### Then call viewMeans() on 'viewslist':
bin_means <- viewMeans(viewslist)
### 'bin_means' is a NumericList object. It as the same "shape" as
### 'mybins' i.e. both are list-like objects of the same length and
### names, and the lengths of 'bin_means[[i]]' and 'mybins[[i]]' are
### the same for any valid i.
### Finally we construct 'averageGranges' from 'mybins' and
'bin_means':
averageGranges <- as(mybins, "GRanges")
mcols(averageGranges)$averagemyValue <- unlist(bin_means)
Hope this helps,
H.
>
> Thanks & kind regards,
> yannick
>
>
>
> ## code is more legible on
>
>
> averageGranges <- GRanges()
> interval_length <- 10000
>
> for (scaffold_index in 1:nrow(myscaffolds)) {
> scaffold <- myscaffolds$id[scaffold_index]
> scaffold_length <- myscaffolds$to[scaffold_index]
>
> num_intervals <- scaffold_length / interval_length
>
> for (interval_index in 0:(num_intervals-1) ) {
> grangesInterval <- GRanges(scaffold, IRanges(interval_index
* interval_length,
> min(
scaffold_length,
>
(interval_index+1)* interval_length )))
> subsetGranges <- subsetByOverlaps( bigGranges ,
grangesInterval)
>
> grangesInterval$averagemyValue <-
mean(subsetGranges$myValue)
> averageGranges <- append(averageGranges, grangesInterval)
> }
> }
>
> # then I use averageGranges
>
> --------------------------------------
> Yannick Wurm
http://yannick.poulet.org
> Ants, Genomes & Evolution ? y.wurm at qmul.ac.uk ? skype:yannickwurm
? +44 207 882 3049
> 5.03A Fogg ? School of Biological & Chemical Sciences ? Queen Mary,
University of London ? Mile End Road ? E1 4NS London ? UK
> Easy custom BLAST interface:
http://www.sequenceserver.com
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at r-project.org
>
https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives:
http://news.gmane.org/gmane.science.biology.informatics.conductor
>
--
Hervé Pagès
Program in Computational Biology
Division of Public Health Sciences
Fred Hutchinson Cancer Research Center
1100 Fairview Ave. N, M1-B514
P.O. Box 19024
Seattle, WA 98109-1024
E-mail: hpages at fhcrc.org
Phone: (206) 667-5791
Fax: (206) 667-1319
