[BioC] A question about the function readGAlignmentPairs in GenomicRnages package

[Niu, Liang (NIH/NIEHS) [E]]

Dear Herve,

Here is the code and output:

> bam.path<-"test.bam"
> bam.file<-BamFile(bam.path,asMates=TRUE)
> 
> param<-ScanBamParam(flag=scanBamFlag(isPaired=TRUE,hasUnmappedMate=FALSE,isUnmappedQuery=FALSE,isNotPassingQualityControls=FALSE,isDuplicate=FALSE,isNotPrimaryRead=FALSE))
> alignment<-readGAlignmentsList(bam.file,param=param)
> length(alignment)
[1] 18041910
> table(mcols(unlist(alignment, use.names=FALSE))$mates[end(PartitioningByEnd(alignment))])

   FALSE 
18041910 

What is the problem?

Thanks!

Liang

________________________________________
From: Hervé Pagès [hpages at fhcrc.org]
Sent: Sunday, March 30, 2014 4:13 PM
To: Niu, Liang (NIH/NIEHS) [E]; bioconductor at r-project.org
Subject: Re: A question about the function readGAlignmentPairs in GenomicRnages package

Dear Liang,

Please keep the communication on the mailing list (by using
the "Reply All" button).

On 03/30/2014 12:13 PM, Niu, Liang (NIH/NIEHS) [E] wrote:
> Dear Herve,
>
> I use the following command to read in the alignments:
>
> bam.path<-"test.bam"
> bam.file<-BamFile(bam.path, yieldSize=1E6, asMates=TRUE)
>
> param<-ScanBamParam(flag=scanBamFlag(isPaired=TRUE,hasUnmappedMate=FALSE,isUnmappedQuery=FALSE,isNotPassingQualityControls=FALSE,isDuplicate=FALSE,isNotPrimaryRead=FALSE))
> alignment<-readGAlignmentsList(bam.file,param=param)
>
> However, each element ((a pair of alignments)) of the GAlignmentsList "alignment"  has meta field Mates as "FALSE" for both reads. Is it a problem?

I'm surprised that *each* list element in your GAlignmentsList object
'alignment' looks like this. How do you know? Why not show us the
object?

My understanding is that even though your BAM file contains paired-end
reads, when you read it with readGAlignmentsList(), you end up with a
GAlignmentsList object where some fraction of the list elements have
the "mates" field set to FALSE. The reads in such list element have the
same QNAME but are not mates. Note that all list elements with the
"mates" field set to TRUE are guaranteed to contain exactly 2 reads.
But there is no such expectation for list elements with the "mates"
field set to FALSE: they can contain 1, 2, 3, or more reads.

So if *each* list element in your GAlignmentsList object
'alignment' really has the "mates" field set to FALSE, that
means none of the reads in your file could be mated. That could
actually happen if your data was not paired-end and if you didn't
use isPaired=TRUE but that doesn't seem to be the case here.

Here is how you can summarize how many list elements have the
"mates" field set to FALSE and how many have it set to TRUE:

   table(mcols(unlist(alignment,
use.names=FALSE))$mates[end(PartitioningByEnd(alignment))])

Finally note that in BioC devel, the "mates" field has been replaced
by the "mate_status" field and that this field is now an attribute of
the list elements rather than of the individual reads. This means that
it's a top-level metadata column (i.e. directly accessible with
mcols(alignment)$mate_status) instead of an inner metadata column
(which are more complicated to access). So it's much easier now to get
the above count:

   table(mcols(alignment)$mate_status)

This will give you something like:

   > table(mcols(galist1)$mate_status)

       mated ambiguous   unmated
       75346         0     21286

As you can see, another change is that this field is now a 3-level
factor (levels are explained in ?readGAlignmentsList).

Hope this helps,
H.


>
> Liang
>
>
> ________________________________________
> From: Hervé Pagès [hpages at fhcrc.org]
> Sent: Sunday, March 30, 2014 2:54 PM
> To: Niu, Liang (NIH/NIEHS) [E]; bioconductor at r-project.org
> Subject: Re: A question about the function readGAlignmentPairs in GenomicRnages package
>
> On 03/30/2014 11:52 AM, Niu, Liang (NIH/NIEHS) [E] wrote:
>> Sorry for the typo in the previous email. What I mean is "I do NEED (not NOT) the pairing information".
>
> I see. So yes, readGAlignmentsList() should do it.
>
> Cheers,
> H.
>
>>
>> Liang
>> ________________________________________
>> From: Hervé Pagès [hpages at fhcrc.org]
>> Sent: Sunday, March 30, 2014 2:51 PM
>> To: Niu, Liang (NIH/NIEHS) [E]; bioconductor at r-project.org
>> Subject: Re: A question about the function readGAlignmentPairs in GenomicRnages package
>>
>> On 03/30/2014 11:31 AM, Niu, Liang (NIH/NIEHS) [E] wrote:
>>> Dear Herve,
>>>
>>> Thanks for your suggestion. I do not the paring information in the downstream analysis, therefore, readGAlignmentsList() is good.
>>
>> Maybe I was not clear enough but if you don't need the pairing, then
>> you can just use readGAlignments(). readGAlignmentsList() will pair
>> everything that can be paired and this has a cost (in terms of
>> performance and memory usage). By using readGAlignments() you avoid
>> that cost and you also end up with an object that is a little bit
>> simpler to manipulate (i.e. a GAlignments instead of GAlignmentsList
>> object).
>>
>> Hope that makes sense,
>> H.
>>
>>>
>>> Liang
>>> ________________________________________
>>> From: Hervé Pagès [hpages at fhcrc.org]
>>> Sent: Sunday, March 30, 2014 2:12 PM
>>> To: Niu, Liang (NIH/NIEHS) [E]
>>> Cc: bioconductor at r-project.org
>>> Subject: Re: A question about the function readGAlignmentPairs in GenomicRnages package
>>>
>>> Hi Liang,
>>>
>>> Hope you don't mind that I'm cc'ing the Bioconductor mailing list, so
>>> other can give suggestions.
>>>
>>> On 03/30/2014 09:32 AM, Niu, Liang (NIH/NIEHS) [E] wrote:
>>>> Dear Mr. Pages,
>>>>
>>>> This is Liang Niu, a research fellow at the National Institute of Environmental Health Sciences.
>>>>
>>>> I am using R to read .bam files for chromatin interaction data sets. Such data sets contains alignments for  paired-end reads from ChIA-PET experiment, thus it has pairs in which two reads on different chromosomes and/or on same strand, and those pairs are valid pairs. I want to use your function readGAlignmentPairs in GenomicRnages package to read  the pairs, but the manual says that it will discard those pairs. Do you have any suggestion?
>>>
>>> You could use readGAlignmentsList() instead of readGAlignmentPairs().
>>> readGAlignmentsList() will keep these "discordant pairs". It will also
>>> keep the reads that cannot be paired.
>>> Note that, depending on what you will do downstream, you don't
>>> necessarily need to pair the reads (e.g. if you're going to compute
>>> the read coverage). In that case you can just load the reads with
>>> readGAlignments().
>>>
>>> Cheers,
>>> H.
>>>
>>>>
>>>> Best,
>>>> Liang
>>>>
>>>
>>> --
>>> 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 at fhcrc.org
>>> Phone:  (206) 667-5791
>>> Fax:    (206) 667-1319
>>>
>>
>> --
>> 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 at fhcrc.org
>> Phone:  (206) 667-5791
>> Fax:    (206) 667-1319
>>
>
> --
> 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 at fhcrc.org
> Phone:  (206) 667-5791
> Fax:    (206) 667-1319
>

--
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 at fhcrc.org
Phone:  (206) 667-5791
Fax:    (206) 667-1319