Michael, Steve,

thanks a lot for your input, this is the functionality i was looking for 
and it works like a charm in the current release, so question solved!!

yet, I have also tried this out in last development version and the 
option 'asRle=TRUE' does not work and gives an error related to IRanges. 
I know the 'devel' version is a moving target but, just in case this was 
unexpected, I've pasted here below the code and its output and the 
session information:

library(BSgenome.Hsapiens.UCSC.hg19)
library(rtracklayer)

download.file("http://hgdownload.cse.ucsc.edu/goldenPath/hg19/phastCons46way/vertebrate/chrY.phastCons46way.wigFix.gz",
               "chrY.phastCons46way.wigFix.gz", method="curl")

si <- Seqinfo(seqnames=seqnames(Hsapiens), seqlengths=seqlengths(Hsapiens))
wigToBigWig("chrY.phastCons46way.wigFix.gz", seqinfo=si)

pc22 <- import.bw(BigWigFile("chrY.phastCons46way.bw")) ## OK!
head(pc22, n=3)
GRanges with 3 ranges and 1 metadata column:
       seqnames         ranges strand |             score
          <Rle>      <IRanges>  <Rle> | <numeric>
   [1]     chrY [10502, 10502]      * |  0.23199999332428
   [2]     chrY [10503, 10503]      * | 0.226999998092651
   [3]     chrY [10504, 10504]      * | 0.218999996781349
   ---
   seqlengths:
        chrY
    59373566

pc22 <- import.bw(BigWigFile("chrY.phastCons46way.bw"), asRle=TRUE) ## 
ERROR!
Error in get(name, envir = asNamespace(pkg), inherits = FALSE) :
   object '.Ranges.coverage' not found

traceback()
12: get(name, envir = asNamespace(pkg), inherits = FALSE)
11: IRanges:::.Ranges.coverage
10: (function (r, v, sl)
     {
         IRanges:::.Ranges.coverage(r, width = sl, weight = v$score)
     })(dots[[1L]][[1L]], dots[[2L]][[1L]], dots[[3L]][[1L]])
9: mapply(function(r, v, sl) {
        IRanges:::.Ranges.coverage(r, width = sl, weight = v$score)
    }, ranges(x), values(x), seqlengths(x), SIMPLIFY = FALSE)
8: .dotargsAsList("Rle", ...)
7: RleList(mapply(function(r, v, sl) {
        IRanges:::.Ranges.coverage(r, width = sl, weight = v$score)
    }, ranges(x), values(x), seqlengths(x), SIMPLIFY = FALSE), compress 
= FALSE)
6: rdToRle(rd)
5: .local(con, format, text, ...)
4: import(con, "BigWig", ...)
3: import(con, "BigWig", ...)
2: import.bw(BigWigFile("chrY.phastCons46way.bw"), asRle = TRUE)
1: import.bw(BigWigFile("chrY.phastCons46way.bw"), asRle = TRUE)

sessionInfo()
R Under development (unstable) (2013-10-20 r64082)
Platform: x86_64-unknown-linux-gnu (64-bit)

locale:
  [1] LC_CTYPE=en_US.UTF8       LC_NUMERIC=C LC_TIME=en_US.UTF8        
LC_COLLATE=en_US.UTF8
  [5] LC_MONETARY=en_US.UTF8    LC_MESSAGES=en_US.UTF8 
LC_PAPER=en_US.UTF8       LC_NAME=C
  [9] LC_ADDRESS=C              LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF8 
LC_IDENTIFICATION=C

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets methods   
base

other attached packages:
  [1] rtracklayer_1.23.1 BSgenome.Hsapiens.UCSC.hg19_1.3.19 BSgenome_1.31.0
  [4] Biostrings_2.31.0 GenomicRanges_1.15.1 XVector_0.3.0
  [7] IRanges_1.21.2 BiocGenerics_0.9.0 BiocInstaller_1.13.1
[10] vimcom_0.9-9 setwidth_1.0-3 colorout_1.0-1

loaded via a namespace (and not attached):
[1] bitops_1.0-6     RCurl_1.95-4.1   Rsamtools_1.15.2 stats4_3.1.0     
tools_3.1.0      XML_3.98-1.1     zlibbioc_1.9.0


thanks again for your help!!
robert.



On 10/23/13 9:43 PM, Michael Lawrence wrote:
>
>
>
> On Wed, Oct 23, 2013 at 9:55 AM, Steve Lianoglou 
> <lianoglou.steve@gene.com <mailto:lianoglou.steve@gene.com>> wrote:
>
>     Hi Robert,
>
>     On Wed, Oct 23, 2013 at 9:03 AM, Robert Castelo
>     <robert.castelo@upf.edu <mailto:robert.castelo@upf.edu>> wrote:
>     > dear list,
>     >
>     > i have to pretty intensively work with genome-wide phastcons
>     scores and
>     > instead of repeatedly interrogate them through the internet via
>     the UCSC
>     > genome browser with 'rtracklayer', i'd prefer to do a bulk
>     download of the
>     > *.phastCons46way.wigFix.gz files (about 0.6Gb) at
>     >
>     >
>     http://hgdownload.cse.ucsc.edu/goldenPath/hg19/phastCons46way/vertebrate
>     >
>     > and then import them into R storing the information in some
>     memory efficient
>     > data structure (Rle?) that provides me also an efficient way to
>     query the
>     > phastcons score at any position of the human genome.
>     >
>     > all documentation and messages i've been able to retrieve
>     through google and
>     > the BioC list correspond to use cases that involve a small
>     fraction of the
>     > genome which can be handled by 'rtracklayer' with an internet
>     connection.
>     >
>     > any hint on how to achieve this goal will be very much appreciated,
>
>     I previously did this by downloading the wig files and converting them
>     to bigWig format. To do so, I used the wigToBigWig tool you can
>     eventually get to form here:
>
>     http://genomewiki.ucsc.edu/index.php/Kent_source_utilities
>
>     It is unfortunate that UCSC doesn't host the bigWig files, as
>     conversion from wig to bigWig is hugely memory intensive (last I
>     remember, anyway).
>
>     (Update: not sure when this happened, but rtracklayer actually wraps
>     this functionality in its wigToBigWig function -- nice! -- so you can
>     convert the wig file to bigWig straight from R (assuming you have a
>     machine with enough horsepower))
>
>
> Of potential importance is that the Kent utilities are licensed for 
> non-commercial use only. That's mostly what drove the implementation 
> of wigToBigWig(). A
>
>     After that, you could then use rtracklayer's import functions over the
>     bigWig files to query them using a GRanges object. Look at the
>     Examples section in the ?import.bw <http://import.bw> man page for
>     some help.
>
>
> In the new release there's an argument asRle that if TRUE will 
> efficiently return an RleList from the BigWig, instead of a GRanges.
>
> If you want to query the RleList for single positions (in a GRanges) 
> and get the result as an atomic vector, it's most efficient to use 
> VariantTools::extractCoverageForPositions. That function should 
> probably live somewhere more general but anyway, it's a fast path for 
> single position RleList=>vector extraction.
>
>     HTH,
>     -steve
>
>     --
>     Steve Lianoglou
>     Computational Biologist
>     Bioinformatics and Computational Biology
>     Genentech
>
>     _______________________________________________
>     Bioconductor mailing list
>     Bioconductor@r-project.org <mailto:Bioconductor@r-project.org>
>     https://stat.ethz.ch/mailman/listinfo/bioconductor
>     Search the archives:
>     http://news.gmane.org/gmane.science.biology.informatics.conductor
>
>


	[[alternative HTML version deleted]]

