[BioC] flip strand of GappedAlignmentPairs Object

Valerie Obenchain vobencha at fhcrc.org
Wed Feb 27 16:48:29 CET 2013


You could use the low-level constructor

     GappedAlignmentPairs(first, last, isProperPair, names=NULL)

'isProperPair' is a logical vector the same length as 'first' and 
'last'.  For example,

     GappedAlignmentPairs(left, right, !logical(length(left)))


Valerie

On 02/27/13 07:12, Stefanie Tauber wrote:
> 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 at fhcrc.org 
> <mailto:vobencha at 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 at r-project.org <mailto:Bioconductor at 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 at univie.ac.at <mailto:stefanie.tauber at univie.ac.at>
> www.cibiv.at <http://www.cibiv.at>
>



More information about the Bioconductor mailing list