[Bioc-sig-seq] GRanges tabulate the coverage

[Patrick Aboyoun]

Dario,
I just checked in a fix to IRanges in BioC 2.6 (IRanges 1.6.13) and BioC 
2.7 (IRanges 1.7.17) that fixes the issue with the table method for 
AtomicList objects. Previous there was just a table,AtomicList-method, 
which was written with CompressedAtomicList objects in mind. I have 
moved the old method definition over to a 
table,CompressedAtomicList-method and have created a new 
table,SimpleAtomicList-method. For your use case, this new method will 
be much faster since it avoids unlisting the data and replicating the 
element names. These changes will become available from bioconductor.org 
within 36 hours.

 > suppressMessages(library(GenomicRanges))
 > library(BSgenome.Hsapiens.UCSC.hg18)
Loading required package: BSgenome
Loading required package: Biostrings

 > chrs <- paste("chr", c(1:22, "X", "Y"), sep = "")
 > myGR <- GRanges(seqnames = c("chr1", "chr3", "chr3", "chr3"), ranges 
= IRanges(100, width = 200), strand = '+')
 > colSums(table(coverage(resize(myGR, 300))))
   0   1   3
198 300 300

 > seqlengths(myGR) <- seqlengths(Hsapiens)[chrs]
 > colSums(table(coverage(resize(myGR, 300))))
          0          1          3
3080418880        300        300

 > sessionInfo()
R version 2.12.0 Under development (unstable) (2010-08-01 r52659)
Platform: i386-apple-darwin9.8.0/i386 (32-bit)

locale:
[1] en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8

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

other attached packages:
[1] BSgenome.Hsapiens.UCSC.hg18_1.3.16 BSgenome_1.17.6
[3] Biostrings_2.17.27                 GenomicRanges_1.1.20
[5] IRanges_1.7.17

loaded via a namespace (and not attached):
[1] Biobase_2.9.0


Cheers,

Patrick


On 8/5/10 11:55 PM, Patrick Aboyoun wrote:
> Dario,
> Thanks for the bug report. This issue is with the 
> table,AtomicList-method when the sum(elementLengths(atomicList)) > 
> 2^31-1. I'll submit a patch to IRanges within the day.
>
> Cheers,
> Patrick
>
>
> On 8/5/10 8:00 PM, Dario Strbenac wrote:
>> Hello,
>>
>> I'm encountering a problem where I can tabulate coverage results when 
>> I've got a GRanges object without seqlengths in it, but when I put 
>> seqlengths in it, it gives an error about negative length vectors. 
>> Here is my small example.
>>
>> library(GenomicRanges)
>> library(BSgenome.Hsapiens.UCSC.hg18)
>>
>> chrs = paste("chr", c(1:22, "X", "Y"), sep = "")
>> myGR<- GRanges(seqnames = c("chr1", "chr3", "chr3", "chr3"), ranges = 
>> IRanges(100, width = 200), strand = '+')
>> coverage<- colSums(table(coverage(resize(myGR, 300))))
>>
>> seqlengths(myGR)<- seqlengths(Hsapiens)[chrs]
>> coverage<- colSums(table(coverage(resize(myGR, 300))))
>>
>> --------------------------------------
>> Dario Strbenac
>> Research Assistant
>> Cancer Epigenetics
>> Garvan Institute of Medical Research
>> Darlinghurst NSW 2010
>> Australia
>>
>> _______________________________________________
>> Bioc-sig-sequencing mailing list
>> Bioc-sig-sequencing at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing
>>
>
> _______________________________________________
> Bioc-sig-sequencing mailing list
> Bioc-sig-sequencing at r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing