[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