I'm working on getting you more information with gdb. In fact I have been processing this data iteratively and in parallel; but sometimes I get lazy and just want coverage() to work ;) While it's true we'll encounter issues above 2 billion, I'm not sure why I would ever need that many reads. 500 million reads is equivalent to 65X coverage WGS for the human genome. Doubling that to 1 billion would be nice for heterogeneous samples; but then returns start diminishing. And reads are getting longer, so we'll need fewer of them. Thanks for dealing with my extreme requests, Michael On Mon, Jun 3, 2013 at 7:48 PM, Martin Morgan <mtmorgan@fhcrc.org> wrote: > On 06/03/2013 07:33 PM, Hervé Pagès wrote: > >> Hi Martin, >> >> On 06/03/2013 06:26 PM, Martin Morgan wrote: >> >>> On 06/03/2013 05:27 PM, Michael Lawrence wrote: >>> >>>> Hey guys, >>>> >>>> Whenever I try to calculate the coverage for a BAM file with more than >>>> say >>>> 500 million reads, I get this error: >>>> >>>> Error in coverage(**readBamGappedAlignments(x, param = param), shift = >>>> shift, : \n error in evaluating the argument 'x' in selecting a method >>>> for function 'coverage': Error in value[[3L]](cond) (from #2) : \n >>>> 'Realloc' could not re-allocate memory (18446744065128005632 bytes)\n >>>> >>>> This looks like integer overflow, possibly within _grow_SCAN_BAM_DATA(). >>>> Could we just use long there? >>>> >>> >>> I wonder if it would be more sensible if less convenient to do this >>> (under Bioc-devel) >>> >>> bf <- open(BamFile(fl, yieldSize=100000000)) >>> cvg <- coverage(**readGAlignmentsFromBam(bf)) >>> while (length(aln <- readGAlignmentsFromBam(bf))) >>> cvg <- cvg + coverage(aln) >>> close(bf) >>> >>> ? It opens the door for better memory management and parallel evaluation. >>> >>> I'm concerned that using size_t (Realloc casts to this) or ptrdiff_t >>> (the size of R long vectors) would only get us through the C code; the >>> representation of this in R would require R long vectors, and Rsamtools >>> does not (yet?) support that. >>> >> >> Sorry if I'm missing something obvious but why would the representation >> of 500 million reads (either as a GappedAlignments object or as a plain >> list as returned by scanBam()) require R long vectors? >> > > not that 500 million would, but that going for 'more' will eventually > (when Michael gets 5 times more ambitious than he is now). > > At the least the software should go up to the limit of R vectors > gracefully; as you point out it shouldn't be having problems with 500 > million reads. > > Martin > > > >> Thanks, >> H. >> >> >> >>> Martin >>> >>> >>>> Michael >>>> >>>> [[alternative HTML version deleted]] >>>> >>>> ______________________________**_________________ >>>> Bioconductor mailing list >>>> Bioconductor@r-project.org >>>> https://stat.ethz.ch/mailman/**listinfo/bioconductor<https: stat="" .ethz.ch="" mailman="" listinfo="" bioconductor=""> >>>> Search the archives: >>>> http://news.gmane.org/gmane.**science.biology.informatics.**condu ctor<http: news.gmane.org="" gmane.science.biology.informatics.conductor=""> >>>> >>>> >>> >>> >> > > -- > Computational Biology / Fred Hutchinson Cancer Research Center > 1100 Fairview Ave. N. > PO Box 19024 Seattle, WA 98109 > > Location: Arnold Building M1 B861 > Phone: (206) 667-2793 > [[alternative HTML version deleted]]