[BioC] [devteam-bioc] readGAlignmentPairsFromBam and readGAlignmentFromBam differences
Emmanouela Repapi
emmanouela.repapi at ludwig.ox.ac.uk
Tue Feb 11 15:55:10 CET 2014
Thank you Herve,
That makes complete sense! Now I get the results I'd expect. I hadn't even seen that ignore.strand parameter because when I did ?findOverlaps I got the IRanges version that doesn't mention the strand at all, and didn't think to look at the GenomicRanges version of it. Thanks a lot!
Emmanouela
Emmanouela Repapi
Computational Biology Research Group
Weatherall Institute of Molecular Medicine
University of Oxford, OX3 9DS
Tel: (01865 2)2253
________________________________________
From: Hervé Pagès [hpages at fhcrc.org]
Sent: 11 February 2014 02:45
To: guest at bioconductor.org; bioconductor at r-project.org; Emmanouela Repapi
Cc: Rsamtools Maintainer
Subject: Re: [devteam-bioc] readGAlignmentPairsFromBam and readGAlignmentFromBam differences
Hi Emmanouela,
Short answer: It's a strand issue (and you should probably use
ignore.strand=TRUE).
Long answer: see below.
On 02/10/2014 07:19 AM, Maintainer wrote:
>
> Hello,
>
> I am doing an RNA-seq analysis for paired-end reads and I have some questions for readGAlignmentPairsFromBam. I have mapped the reads with tophatv2 and aligned them with Genomic Features using the readGAlignmentPairsFromBam function, see the code below.
> I understand that since I have paired-end reads I should use readGAlignmentPairsFromBam instead of readGAlignmentFromBam but it seems to have much more strict criteria for the pairing and therefore, the counts are much lower (mean number of counts: 185 versus 440 for a random sample).
>
> I read the pairing criteria trying to understand why lots of the reads are not being counted but I get differences for which the reason of exclusion is not clear to me.
>
> example for one of these reads:
> $ samtools view 209_accepted_hits.bam | grep "HISEQ2000-02:440:C35VLACXX:7:1102:3519:68259"
> HISEQ2000-02:440:C35VLACXX:7:1102:3519:68259 97 chr15 77414942 50 51M = 77415180 289 GTCTGCAGAGTCTCCACTCTGCCGGAAGTCGGATCCCCTTTTAGACTCCAG &&&(((((**(**++++++++++++++++++++)+++++++++++++++++ AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:51 YT:Z:UU XS:A:- NH:i:1
> HISEQ2000-02:440:C35VLACXX:7:1102:3519:68259 145 chr15 77415180 50 51M = 77414942 -289 TGCATTTTGATCCAGTCACCGTTTCTTTAGGCTGCTCTGCCCTTAACCCAG *+(#++)+)+++**($)(++++*)(+++*)'++++*)+**(((&&%%&$$$ AS:i:-3 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:20A30 YT:Z:UU XS:A:- NH:i:1
If I put those 2 records alone in a BAM file, they show up with
readGAlignmentsFromBam():
> gal1 <- readGAlignmentsFromBam("test1.bam", use.names=TRUE)
> gal1
GAlignments with 2 alignments and 0 metadata columns:
seqnames strand
cigar qwidth start end
<Rle> <Rle>
<character> <integer> <integer> <integer>
HISEQ2000-02:440:C35VLACXX:7:1102:3519:68259 chr15 +
51M 51 77414942 77414992
HISEQ2000-02:440:C35VLACXX:7:1102:3519:68259 chr15 -
51M 51 77415180 77415230
width ngap
<integer> <integer>
HISEQ2000-02:440:C35VLACXX:7:1102:3519:68259 51 0
HISEQ2000-02:440:C35VLACXX:7:1102:3519:68259 51 0
---
seqlengths:
chr15
103494974
They also get paired by readGAlignmentPairsFromBam():
> galp1 <- readGAlignmentPairsFromBam("test1.bam", use.names=TRUE)
> galp1
GAlignmentPairs with 1 alignment pair and 0 metadata columns:
seqnames strand :
ranges -- ranges
<Rle> <Rle> :
<IRanges> -- <IRanges>
HISEQ2000-02:440:C35VLACXX:7:1102:3519:68259 chr15 + :
[77414942, 77414992] -- [77415180, 77415230]
---
seqlengths:
chr15
103494974
>
> is included in the result from readGAlignmentFromBam but not from readGAlignmentPairsFromBam. Also, many pairs with bitwise flags of 97 and 145 (for pair) are not included in readGAlignmentPairsFromBam although their flags satisfy the pairing criteria.
>
> Also I get some reads in readGAlignmentPairsFromBam that I don’t get from readGAlignmentFromBam, which again I don’t understand why.
>
> Example for one of these reads:
> $samtools view 209_accepted_hits.bam.fltrd.bam | grep "HISEQ2000-02:440:C35VLACXX:7:2311:15754:64830"
> HISEQ2000-02:440:C35VLACXX:7:2311:15754:64830 161 chr15 77426477 50 51M = 77426641 215 TGTGGTTCCAGTCACCAGCTGCACATCTGAGACACACCAGACACAAGCTCT %$%&('&(**)(*+++)++++++++++++++++++++++++++++++++++ AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:51 YT:Z:UU XS:A:- NH:i:1
> HISEQ2000-02:440:C35VLACXX:7:2311:15754:64830 81 chr15 77426641 50 51M = 77426477 -215 GCCATAGCTGTCCATGATCTCTGATTACCTTTCTTCTATGATTTAAAAGGG ++++++++++++++++++++++++++*+++++++++++*****(((((&&& AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:51 YT:Z:UU NH:i:1
Here too, if I put those 2 records alone in a BAM file:
> gal2 <- readGAlignmentsFromBam("test2.bam", use.names=TRUE)
> gal2
GAlignments with 2 alignments and 0 metadata columns:
seqnames strand
cigar qwidth start end
<Rle> <Rle>
<character> <integer> <integer> <integer>
HISEQ2000-02:440:C35VLACXX:7:2311:15754:64830 chr15 +
51M 51 77426477 77426527
HISEQ2000-02:440:C35VLACXX:7:2311:15754:64830 chr15 -
51M 51 77426641 77426691
width ngap
<integer> <integer>
HISEQ2000-02:440:C35VLACXX:7:2311:15754:64830 51 0
HISEQ2000-02:440:C35VLACXX:7:2311:15754:64830 51 0
---
seqlengths:
chr15
103494974
Also:
> galp2 <- readGAlignmentPairsFromBam("test2.bam", use.names=TRUE)
> galp2
GAlignmentPairs with 1 alignment pair and 0 metadata columns:
seqnames strand :
ranges -- ranges
<Rle> <Rle> :
<IRanges> -- <IRanges>
HISEQ2000-02:440:C35VLACXX:7:2311:15754:64830 chr15 - :
[77426641, 77426691] -- [77426477, 77426527]
---
seqlengths:
chr15
103494974
So the "count" issue you describe is not a "pairing criteria" issue.
You need to understand how findOverlaps() works on a GAlignmentPairs
object. First with GAlignments object 'gal1':
> findOverlaps(range.exon["328561"], gal1)
Hits of length 1
queryLength: 1
subjectLength: 2
queryHits subjectHits
<integer> <integer>
1 1 2
Your gene (which is on the minus strand) overlaps with the 2nd element
in 'gal1'.
When you pass a GAlignments or GAlignmentPairs object to findOverlaps(),
it's first turned into a GRangesList object internally, with grglist().
The key here is that, grglist() on a GAlignmentPairs object inverts the
strand of the 2nd mate:
> grglist(galp1)
GRangesList of length 1:
$HISEQ2000-02:440:C35VLACXX:7:1102:3519:68259
GRanges with 2 ranges and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
[1] chr15 [77414942, 77414992] +
[2] chr15 [77415180, 77415230] +
---
seqlengths:
chr15
103494974
See ?GAlignmentPairs for more info about this (especially the IMPORTANT
note in the Coercion section). The reason for this behavior is because
the strand of the 2nd mate was already inverted by the sequencing
technology (when the cDNA template was sequenced from the 2nd end).
This is why you don't see any overlap when you do:
> length(findOverlaps(range.exon["328561"], galp1))
[1] 0
(because your gene is on the minus strand).
Most of the time those strand considerations don't matter because the
RNA-seq protocol is not stranded. In that case you should use
ignore.strand=TRUE:
> findOverlaps(range.exon["328561"], gal1, ignore.strand=TRUE)
Hits of length 2
queryLength: 1
subjectLength: 2
queryHits subjectHits
<integer> <integer>
1 1 1
2 1 2
> findOverlaps(range.exon["328561"], galp1, ignore.strand=TRUE)
Hits of length 1
queryLength: 1
subjectLength: 1
queryHits subjectHits
<integer> <integer>
1 1 1
Similar story for gal2/galp2:
> length(findOverlaps(range.exon["328561"], gal2))
[1] 0
No hit with the 1st element in 'gal2' because even though gal2[1]
is in the range of your gene, it's not on the same strand.
No hit with the 2nd element in 'gal2' because even though gal2[2]
is on the same strand as your gene, it's not in its range.
However, we see a hit with the pair:
> length(findOverlaps(range.exon["328561"], galp2))
[1] 1
That's because when converting to GRangesList, everything ends up
on the minus strand:
> grglist(galp2)
GRangesList of length 1:
$HISEQ2000-02:440:C35VLACXX:7:2311:15754:64830
GRanges with 2 ranges and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
[1] chr15 [77426477, 77426527] -
[2] chr15 [77426641, 77426691] -
---
seqlengths:
chr15
103494974
But again, if you use ignore.strand=TRUE:
> findOverlaps(range.exon["328561"], gal2, ignore.strand=TRUE)
Hits of length 1
queryLength: 1
subjectLength: 2
queryHits subjectHits
<integer> <integer>
1 1 1
Hope this helps,
H.
>
> Can someone, please, explain why this is?
> Thank you in advance for your time.
> Emmanouela
>
>
> -- output of sessionInfo():
>
>
>
> library(GenomicRanges)
> library(Rsamtools)
> library(GenomicFeatures)
> library(rtracklayer)
>
> txdb.ensgene <- makeTranscriptDbFromUCSC(genome = "mm9", tablename = "refGene")
> range.exon <- exonsBy(txdb.ensgene, "gene")
>
> my_bam <- "209_accepted_hits.bam.fltrd.bam"
> aligns <- readGAlignmentPairsFromBam(my_bam,use.names=T)
> aligns_s <- readGAlignmentsFromBam(my_bam,use.names=T)
>
> temp3 <- findOverlaps(range.exon["328561"],aligns) # gene Apol10b
> temp3_aligns <- aligns[subjectHits(temp3)]
>
> temp4 <- findOverlaps(range.exon["328561"],aligns_s)
> temp4_aligns <- aligns_s[subjectHits(temp4)]
>
> setdiff(names(temp3_aligns) , names(temp4_aligns))[1]
> [1] "HISEQ2000-02:440:C35VLACXX:7:2311:15754:64830"
>
> setdiff(names(temp4_aligns),names(temp3_aligns))[1]
> [1] "HISEQ2000-02:440:C35VLACXX:7:1102:3519:68259"
>
>
>
> sessionInfo()
> R version 3.0.1 (2013-05-16)
> Platform: x86_64-unknown-linux-gnu (64-bit)
>
> locale:
> [1] LC_CTYPE=en_GB.ISO-8859-1 LC_NUMERIC=C
> [3] LC_TIME=en_GB.ISO-8859-1 LC_COLLATE=en_GB.ISO-8859-1
> [5] LC_MONETARY=en_GB.ISO-8859-1 LC_MESSAGES=en_GB.ISO-8859-1
> [7] LC_PAPER=C LC_NAME=C
> [9] LC_ADDRESS=C LC_TELEPHONE=C
> [11] LC_MEASUREMENT=en_GB.ISO-8859-1 LC_IDENTIFICATION=C
>
> attached base packages:
> [1] parallel stats graphics grDevices utils datasets methods
> [8] base
>
> other attached packages:
> [1] org.Mm.eg.db_2.10.1 RSQLite_0.11.4 DBI_0.2-7
> [4] annotate_1.40.0 rtracklayer_1.22.1 GenomicFeatures_1.14.2
> [7] AnnotationDbi_1.24.0 Biobase_2.22.0 Rsamtools_1.14.2
> [10] Biostrings_2.30.1 GenomicRanges_1.14.4 XVector_0.2.0
> [13] IRanges_1.20.6 BiocGenerics_0.8.0
>
> loaded via a namespace (and not attached):
> [1] biomaRt_2.18.0 bitops_1.0-6 BSgenome_1.30.0 RCurl_1.95-4.1
> [5] stats4_3.0.1 tools_3.0.1 XML_3.98-1.1 xtable_1.7-1
> [9] zlibbioc_1.8.0
>
>
> --
> Sent via the guest posting facility at bioconductor.org.
>
>
>
> ________________________________________________________________________
> devteam-bioc mailing list
> To unsubscribe from this mailing list send a blank email to
> devteam-bioc-leave at lists.fhcrc.org
> You can also unsubscribe or change your personal options at
> https://lists.fhcrc.org/mailman/listinfo/devteam-bioc
>
--
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
More information about the Bioconductor
mailing list