I have observed big slowdowns with GRanges when the number of seqlevels is
large; not sure if this is related.  My observation was mostly about
findOverlaps.


On Wed, May 8, 2013 at 11:27 AM, Stefanie Tauber <
stefanie.tauber@univie.ac.at> wrote:

> To be more specific,
> it takes about 10 minutes to calculate the coverage of 1 x 10^6 reads on
> my pc (linux, 4 cpus, 8GB RAM)
>
>
> This seems a bit long to me, what do you think?
>
>
>
> Am 08.05.2013 um 17:09 schrieb Steve Lianoglou <lianoglou.steve@gene.com>:
>
> > Another approach to what Martin suggests would be to simply read in the
> reads and process them (do the coverage calculation per transcript stuff)
> one chromosome at a time.
> >
> > On Wednesday, May 8, 2013, Martin Morgan wrote:
> > On 05/08/2013 04:40 AM, Stefanie Tauber wrote:
> > Dear list,
> >
> > I have some human data, mapped with bowtie to the transcriptome and then
> > further processed with eXpress.
> > The output of eXpress is a sorted bam file, containing the alignments of
> > the reads to different transcripts.
> >
> > So, the bam file contains about 48 x 10^6 alignments, distributed over
> > about 50000 transcripts.
> > I read in the bam file with 'readGappedAlignments'.
> >
> >
> > reads
> > GappedAlignments with 46765032 alignments and 0 metadata columns:
> >                                   seqnames strand       cigar    qwidth
> >                                      <Rle>  <Rle> <character> <integer>
> >           [1] hg19_ensGene_ENST00000237247      +         75M        75
> >           [2] hg19_ensGene_ENST00000371036      -         75M        75
> >           [3] hg19_ensGene_ENST00000371037      -         75M        75
> >           [4] hg19_ensGene_ENST00000471889      -         75M        75
> >           [5] hg19_ensGene_ENST00000377479      +         75M        75
> >           ...                          ...    ...         ...       ...
> >    [46765028] hg19_ensGene_ENST00000344867      -         75M        75
> >    [46765029] hg19_ensGene_ENST00000400848      -         75M        75
> >    [46765030] hg19_ensGene_ENST00000400848      -         75M        75
> >    [46765031] hg19_ensGene_ENST00000400848      -         75M        75
> >    [46765032] hg19_ensGene_ENST00000400829      -         75M        75
> >                   start       end     width      ngap
> >               <integer> <integer> <integer> <integer>
> >           [1]      1783      1857        75         0
> >           [2]       690       764        75         0
> >           [3]      1632      1706        75         0
> >           [4]      2121      2195        75         0
> >           [5]       391       465        75         0
> >           ...       ...       ...       ...       ...
> >    [46765028]       280       354        75         0
> >    [46765029]       509       583        75         0
> >    [46765030]       509       583        75         0
> >    [46765031]       509       583        75         0
> >    [46765032]       179       253        75         0
> >    ---
> >    seqlengths:
> >     hg19_ensGene_ENST00000237247 ... hg19_ensGene_ENST00000400829
> >                             2580 ...                          999
> >
> >
> >
> >
> > Now, for whatever reason, I would like to have the coverage vector for
> > each transcript.
> > So, basically I would just compute
> > coverage(reads).
> >
> > This works very well for smaller genomes and experiments (e.g. yeast with
> > 10 x 10^6 alignments).
> >
> > But, using 'coverage' on this large ReadGappedAlignments Object takes
> > really very very long.
> >
> > Hi Stefanie --
> >
> > One reason might be that you are on the edge of physical memory
> availability and you are swapping to disk. A solution is to read in chunks,
> which turns out to be very easy
> >
> >   bam = BamFile(<path/to/bam>, yieldSize = 10000000)
> >   open(bam)
> >   cvg = NULL
> >   while (length(reads <- readGappedAlignments(bam, <...>))) {
> >       if (cvg == NULL)
> >           cvg = coverage(reads)
> >       else cvg = cvg + coverage(reads)
> >   }
> >   close(bam)
> >
> > This will be relatively efficient; yieldSize can be adjusted, including
> to a size so that several bam files can be processed in parallel. Start
> with a small yieldSize to quickly debug the iteration, before processing
> the whole file. It can be further simplified using coverage,BamFile-method
> directly, if you're not wanting to do anything else with the gapped
> alignments.
> >
> > Martin
> >
> >
> > Is there a way to do this better (in R)? Or is this just not feasible to
> > do in R and I should rather preprocess the data outside R?
> >
> > I would appreciate any comments,
> > best,
> > Stefanie
> >
> > _______________________________________________
> > Bioconductor mailing list
> > Bioconductor@r-project.org
> > https://stat.ethz.ch/mailman/listinfo/bioconductor
> > Search the archives:
> 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
> >
> > _______________________________________________
> > Bioconductor mailing list
> > Bioconductor@r-project.org
> > https://stat.ethz.ch/mailman/listinfo/bioconductor
> > Search the archives:
> http://news.gmane.org/gmane.science.biology.informatics.conductor
> >
> >
> > --
> > Steve Lianoglou
> > Computational Biologist
> > Department of Bioinformatics and Computational Biology
> > Genentech
> >
>
> DI Stefanie Tauber
>
> Center for Integrative Bioinformatics Vienna (CIBIV)
> (CIBIV is a joint institute of Vienna University and Medical University)
> Max F. Perutz Laboratories (MFPL)
> Campus Vienna Biocenter 5 (VBC5), Ebene 1, Room 1812.2
> Dr. Bohr Gasse 9
> A-1030 Wien, Austria
> Phone: ++43 +1 / 42772-4030
> Fax:     ++43 +1 / 42772-4098
> email:   stefanie.tauber@univie.ac.at
> www.cibiv.at
>
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> Bioconductor mailing list
> 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]]

