[BioC] flip strand of GappedAlignmentPairs Object
Valerie Obenchain
vobencha at fhcrc.org
Wed Feb 27 16:06:17 CET 2013
Hi Stephanie,
A GappedAlignmentPairs object is made up of a 'left' and a 'right'
GappedAlignments.
> ex1_file <- system.file("extdata", "ex1.bam", package="Rsamtools")
> galp <- readGappedAlignmentPairs(ex1_file, use.names=TRUE)
You can extract the left/right parts with an accessor then look at
strand, cigars, gaps etc.
> left <- left(galp)
> class(left)
[1] "GappedAlignments"
attr(,"package")
[1] "GenomicRanges"
> strand(left)
factor-Rle of length 1572 with 665 runs
Lengths: 3 25 1 1 1 3 1 1 3 2 2 ... 1 6 1 2 3 2 1 5
1 51
Values : - + - + - + - + - + - ... + - + - + - + -
+ -
Levels(3): + - *
> head(cigar(left))
[1] "35M" "36M" "35M" "36M" "35M" "35M"
> head(ngap(left))
[1] 0 0 0 0 0 0
It's on these left/right pairs that you can reset the strand.
> strand(left)[1:3]
factor-Rle of length 3 with 1 run
Lengths: 3
Values : +
Levels(3): + - *
> strand(left)[1:3] <- "-"
> strand(left)[1:3]
factor-Rle of length 3 with 1 run
Lengths: 3
Values : -
Levels(3): + - *
Please be sure to read the details on the man page regarding how the
strands were assigned.
?GappedAlignmentPairs
Valerie
On 02/27/13 05:55, Stefanie Tauber wrote:
> Dear list,
>
> Let's say I have single-end strand-specific RNA-Seq data.
> Then I could read in the data with:
>
> library(GenomicRanges)
>
>> reads = readGappedAlignments("data.bam")
> As I typically don't know which strand has been sequenced, it might be necessary to flip the strand of the reads before doing any kind
> of summarizeOverlap:
>
>> strand(reads) = ifelse(strand(reads) == "+", "-", "+")
> Now I could do a 'summarizeOverlaps' between the reads and any kind of transcript database.
>
> Now, if I have paired-end strand-specific RNA-Seq data.
> I read the data:
>> reads = readGappedAlignmentPairs("data.bam")
> How can I flip the strand of the reads? As it's now paired-end data, it does not work as above.
> I could imagine several work-arounds, I just wonder if there is any straight approach which I miss at the moment..
>
>
> best regards,
> Stefanie
>
>
>
>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor
More information about the Bioconductor
mailing list