Hi Valerie,

thanks for the quick answer.
I  have already checked how the strand is assigned for paired-end data.
After setting the strand for the left and right parts, 
can I somehow rebuild the gappedalignmentPairs object?

best,
Stefanie

Am 27.02.2013 um 16:06 schrieb Valerie Obenchain <vobencha@fhcrc.org>:

> 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@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/bioconductor
>> Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor
> 

DI Stefanie Tauber

Center for Integrative Bioinformatics Vienna (CIBIV)
(CIBIV is a joint institute of Vienna University and Medical University)
Max F. Perutz Laboratories (MFPL)
Campus Vienna Biocenter 5 (VBC5), Ebene 1, Room 1812.2
Dr. Bohr Gasse 9
A-1030 Wien, Austria
Phone: ++43 +1 / 42772-4030
Fax:     ++43 +1 / 42772-4098
email:   stefanie.tauber@univie.ac.at
www.cibiv.at


	[[alternative HTML version deleted]]

