[BioC] flip strand information in GappedAlignments or GRangesList Object
Nicolas Delhomme
delhomme at embl.de
Fri Mar 30 21:04:04 CEST 2012
Hi Hervé,
On 30 Mar 2012, at 20:57, Hervé Pagès wrote:
> Hi Stefanie,
>
> On 03/30/2012 07:09 AM, Michael Lawrence wrote:
>> On Fri, Mar 30, 2012 at 5:54 AM, Steve Lianoglou<
>> mailinglist.honeypot at gmail.com> wrote:
>>
>>> Hi,
>>>
>>> On Fri, Mar 30, 2012 at 8:42 AM, Kathi Zarnack<zarnack at ebi.ac.uk> wrote:
>>>> Hi,
>>>> I think you have to add vector() to strand(ga), since the strand
>>> information
>>>> comes as Rle which does not work in the ifelse() test directly. At least
>>>> this is my experience with GRanges objects.
>>>
>>> It's true that sometimes you get bitten by Rle's popping up where you
>>> didn't expect them when you try and index things and need to do some
>>> as.vector() mojo here and there, but in this case it should actually
>>> work (depending on your version of R/bioc, I guess).
>>>
>>> In my case (running R-2.15 RC), the Rle "surprise" doesn't catch you
>>> here, but it's still a good idea to keep in mind.
>>>
>>>
>> Thanks, Steve. ifelse() has worked on Rle for a very long time. Btw, this
>> issue usually comes up when trying to extract from an ordinary vector with
>> [. Since we can't / don't want to define methods on base generics for base
>> data structures, we need to define new generics. In the case of [, one can
>> use seqselect() as an alternative.
>>
>> Thanks for pointing that out!
>
> And in case you wonder how to do this on a GRangesList object, the
> idiom is:
>
> (1) Unlist:
>
> unlisted <- unlist(grl, use.names=FALSE)
>
> (2) Applies Steve's ifelse() code to 'unlisted' (GRanges object).
>
> (3) Relist:
>
> grl2 <- relist(unlisted, grl)
>
> The "unlist/transform/relist" idiom is a very efficient/convenient way
> to transform CompressedList objects in general (not only GRangesList
> objects, which are only a particular case of CompressedList objects).
> However, it does NOT work for any transformation performed in (2),
> but only for a transformation where:
> (a) the input/output have the same class
> (b) and they have the same length (and the i-th element in the
> output corresponds to the i-th element in the input, e.g. rev()
> does not satisfy this).
>
> Note that relist() looses the original elementMetadata() but you can
> always propagate it "by hand":
>
> elementMetadata(grl2) <- elementMetadata(grl)
>
> Maybe we should modify the "relist" methods to do this automatically,
> after all it propagates the names so why not the elementMetadata...
>
I did not know about relist, so I was doing it my own way (which I had trouble making). It would indeed be excellent if relist would propagate the metadata as well!
Cheers,
N.
> Cheers,
> H.
>
>>>
>>> -steve
>>>
>>>> Best regards,
>>>> Kathi
>>>>
>>>>
>>>> On 30/03/12 13:37, Steve Lianoglou wrote:
>>>>>
>>>>> Hi,
>>>>>
>>>>> On Fri, Mar 30, 2012 at 6:04 AM, Stefanie<stefanie.tauber at univie.ac.at>
>>>>> wrote:
>>>>>>
>>>>>> Dear list,
>>>>>>
>>>>>> having a GappedAlignments or GRangesList object at hand,
>>>>>> what is the quickest way to flip strand signs?
>>>>>>
>>>>>> So for each entry, I want to flip "-" to "+" and vice versa.
>>>>>
>>>>> Let's assume that your GappedAlignments object is called `ga`, I think
>>>>> this should work, no?
>>>>>
>>>>> R> strand(ga)<- ifelse(strand(ga) == '+', '-', '+')
>>>>>
>>>>> To the devs:
>>>>>
>>>>> For some reason, `example(GappedAlignments)` is throwing the following
>>>>> error on me, so I can't actually test at the moment:
>>>>>
>>>>> Error in elementLengths(rglist(x)) :
>>>>> error in evaluating the argument 'x' in selecting a method for
>>>>> function 'elementLengths': Error in .Call(.NAME, ..., PACKAGE =
>>>>> PACKAGE) :
>>>>> Incorrect number of arguments (6), expecting 4 for
>>>>> 'cigar_to_list_of_IRanges_by_alignment'
>>>>>
>>>>> sessionInfo() pasted below.
>>>>>
>>>>> HTH,
>>>>>
>>>>> -steve
>>>>>
>>>>> R version 2.15.0 RC (2012-03-24 r58823)
>>>>> Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)
>>>>>
>>>>> locale:
>>>>> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
>>>>>
>>>>> attached base packages:
>>>>> [1] stats graphics grDevices utils datasets methods base
>>>>>
>>>>> other attached packages:
>>>>> [1] Rsamtools_1.7.41 Biostrings_2.23.6 GenomicRanges_1.7.40
>>>>> [4] IRanges_1.13.34 BiocGenerics_0.1.14
>>>>>
>>>>
>>>> --
>>>> Dr. Kathi Zarnack
>>>> Luscombe Group
>>>> European Bioinformatics Institute
>>>> Wellcome Trust Genome Campus
>>>> Hinxton, Cambridge
>>>> CB10 1SD, UK
>>>> tel +44 1223 494 526
>>>>
>>>
>>>
>>>
>>> --
>>> Steve Lianoglou
>>> Graduate Student: Computational Systems Biology
>>> | Memorial Sloan-Kettering Cancer Center
>>> | Weill Medical College of Cornell University
>>> Contact Info: http://cbio.mskcc.org/~lianos/contact
>>>
>>> _______________________________________________
>>> 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
>>>
>>
>> [[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
>
>
> --
> Hervé Pagès
>
> Program in Computational Biology
> Division of Public Health Sciences
> Fred Hutchinson Cancer Research Center
> 1100 Fairview Ave. N, M1-B514
> P.O. Box 19024
> Seattle, WA 98109-1024
>
> E-mail: hpages at fhcrc.org
> Phone: (206) 667-5791
> Fax: (206) 667-1319
>
> _______________________________________________
> 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