Hi Herve,
Thanks! I have a server with 96GB of memory, so now swap issues.
I just downloaded and installed the dev branch of bioconductor on my server
so that I have the latest Rsamtools to work with.

We are using "Star" to map the reads and it looks like it produces many
multi-mapped reads and the mapping is to > 2 reads. I could see that
"readGappedAlignmentPairs"  was working very hard.... a lot of diagnostic
messges about mates being dropped.

Thanks for the tip on the coverage using grlist.
-best
-Lax
On Fri, Jun 22, 2012 at 2:11 PM, Hervé Pagès <hpages@fhcrc.org> wrote:

> Hi Iyer,
>
>
> On 06/22/2012 06:56 AM, Lakshmanan Iyer wrote:
>
>> Hi
>> My apologies for multiple posting if it happens-- I sent the last mails
>> from other accounts which may not be registered with Bioc-list
>> -Lax
>> Two questions:
>>
>> 1. Is readGapppedAlignmentPairs - the most efficient way to read a
>> paired-end bam file with mulit-mapped reads?
>> I am asking as it takes an enormous amount of time to process and load.
>>
>
> With a recent version of Rsamtools (>= 1.9.10), last time I tried it
> took between 30 and 40 min to call readGapppedAlignmentPairs() on a BAM
> file with 70 million records (35 million pairs). Even for a fixed nb of
> records and using the same machine, timing could vary a lot depending
> on how "hard" it is to pair the records which mostly depends on the
> average number of records sharing the same QNAME (the lowest, the
> easiest). You can get that number using the new quickBamCounts() utility
> from Rsamtools:
>
> > quickBamCounts("AML_330-0_**gsnap_filter.bam")
>                                group |    nb of |    nb of | mean / max
>                                   of |  records |   unique | records per
>                              records | in group |   QNAMEs | unique QNAME
> All records.......................**. A | 70446309 | 35407913 | 1.99 / 2
>  o template has single segment.... S |        0 |        0 |   NA / NA
>  o template has multiple segments. M | 70446309 | 35407913 | 1.99 / 2
>      - first segment.............. F | 35313360 | 35313360 |    1 / 1
>      - last segment............... L | 35132949 | 35132949 |    1 / 1
>      - other segment.............. O |        0 |        0 |   NA / NA
>
> Note that (S, M) is a partitioning of A, and (F, L, O) is a partitioning
> of M.
> Indentation reflects this.
>
> Here the average number of records per unique QNAME is 1.99 (and the max
> is 2), which is ideal.
>
> I don't remember the exact amount of memory needed to load that file
> with 70M records though. All I remember is that you need a lot :-/
> (probably > 20GB). Make sure you have enough memory and that your system
> is not swapping.
>
> readGapppedAlignmentPairs() is basically calling readGapppedAlignments()
> followed by findMateAlignment(). Each of the 2 operations are expensive
> and IIRC I'm not sure there are obvious/easy optimizations that could
> be done on readGapppedAlignments() itself. However, for
> findMateAlignment(), there are a few easy ones that have been
> on my list for a couple of months now and that I will implement ASAP.
> They should improve speed and reduce memory usage.
>
>
>
>> 2. How does one work with coverage on GappedAlignmentPairs  in the context
>> of RNASeq?
>> The simplest way is to consider each left and right read as separate -
>> essentially loose the "paired" information and calculate coverage.
>> However, if both the left and right pair reads fall within a feature of
>> interest - say an exon, does it imply coverage of the region of the exon
>> between the reads too
>>
>>
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>>>>>>>>>>>    LLLLLLLLLL--------------------**
>> -------RRRRRRRRRR
>>                     ^^^^^^^^^^^^^^^^^
>>
>> In the figure above, the exon is represented by ">" and L and R represents
>> the left and right reads aligned to the exon.
>> I am talking about the region represented by "^".  Do we assume coverage
>> for this region too?
>> Does Coverage on GappedAlignmentPairs do this?
>>
>
> No, coverage on a GappedAlignmentPairs does not do this, but
> 'coverage(range(grglist(galp))**)' will do this.
>
> Cheers,
> H.
>
>
>> -best
>> -Lax
>> Center for Neuroscience Research
>> Tufts Univeristy School of Medicine
>> Boston, MA
>>
>>        [[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.**conductor<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@fhcrc.org
> Phone:  (206) 667-5791
> Fax:    (206) 667-1319
>

	[[alternative HTML version deleted]]

