[BioC] flip strand of GappedAlignmentPairs Object

[Valerie Obenchain]

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