[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