flip strand information in GappedAlignments or GRangesList Object
1
1
Entering edit mode
Stefanie ▴ 360
@stefanie-5192
Last seen 10.2 years ago
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. Thanks, Stefanie
• 1.7k views
ADD COMMENT •
1
Entering edit mode
@steve-lianoglou-2771
Last seen 20 months ago
United States
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 -- 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
ADD COMMENT •
0
Entering edit mode
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. 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
ADD REPLY •
0
Entering edit mode
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 for pointing that out! -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
ADD REPLY •
0
Entering edit mode
On Fri, Mar 30, 2012 at 5:54 AM, Steve Lianoglou < mailinglist.honeypot@gmail.com> wrote: > Hi, > > On Fri, Mar 30, 2012 at 8:42 AM, Kathi Zarnack <zarnack@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! > > -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@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@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]]
ADD REPLY •
0
Entering edit mode
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... 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
ADD REPLY •
0
Entering edit mode
One note about propagating metadata, I sometimes use relist to form a list out of an arbitrary flat structure that should have the same structure as an existing list. In this use, I am not sure it would make sense to carry over the metadata. Maybe as an option? Michael On Fri, Mar 30, 2012 at 11:57 AM, Hervé Pagès <hpages@fhcrc.org> 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@gmail.com**> wrote: >> >> Hi, >>> >>> On Fri, Mar 30, 2012 at 8:42 AM, Kathi Zarnack<zarnack@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... > > 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@**>>>>> univie.ac.at <stefanie.tauber@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<http: cbio.="" mskcc.org="" %7elianos="" contact=""> >>> >>> ______________________________**_________________ >>> Bioconductor mailing list >>> Bioconductor@r-project.org >>> https://stat.ethz.ch/mailman/**listinfo/bioconductor<https: stat.="" ethz.ch="" mailman="" listinfo="" bioconductor=""> >>> Search the archives: >>> http://news.gmane.org/gmane.**science.biology.informatics.**conduc tor<http: news.gmane.org="" gmane.science.biology.informatics.conductor=""> >>> >>> >> [[alternative HTML version deleted]] >> >> >> ______________________________**_________________ >> Bioconductor mailing list >> Bioconductor@r-project.org >> https://stat.ethz.ch/mailman/**listinfo/bioconductor<https: stat.e="" thz.ch="" mailman="" listinfo="" bioconductor=""> >> Search the archives: http://news.gmane.org/gmane.** >> science.biology.informatics.**conductor<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@fhcrc.org > Phone: (206) 667-5791 > Fax: (206) 667-1319 > [[alternative HTML version deleted]]
ADD REPLY •
0
Entering edit mode
Hi Michael, On 03/30/2012 02:56 PM, Michael Lawrence wrote: > One note about propagating metadata, I sometimes use relist to form a > list out of an arbitrary flat structure that should have the same > structure as an existing list. In this use, I am not sure it would make > sense to carry over the metadata. Maybe as an option? I don't see much difference with the names of the skeleton, which are carried over. More generally speaking I think we should try to treat names and elementMetdata more consistently. H. > > Michael > > On Fri, Mar 30, 2012 at 11:57 AM, Hervé Pagès <hpages at="" fhcrc.org=""> <mailto:hpages at="" fhcrc.org="">> 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 > <mailto:mailinglist.honeypot at="" gmail.com="">__> wrote: > > Hi, > > On Fri, Mar 30, 2012 at 8:42 AM, Kathi > Zarnack<zarnack at="" ebi.ac.uk="" <mailto: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... > > 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=""> <mailto: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 <tel:%2b44%201223%20494%20526> > > > > > -- > 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 > <http: cbio.mskcc.org="" %7elianos="" contact=""> > > _________________________________________________ > Bioconductor mailing list > Bioconductor at r-project.org <mailto:bioconductor at="" r-project.org=""> > https://stat.ethz.ch/mailman/__listinfo/bioconductor > <https: stat.ethz.ch="" mailman="" listinfo="" bioconductor=""> > Search the archives: > http://news.gmane.org/gmane.__science.biology.informatics.__conductor > <http: news.gmane.org="" gmane.science.biology.informatics.conductor=""> > > > [[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 > <https: stat.ethz.ch="" mailman="" listinfo="" bioconductor=""> > Search the archives: > http://news.gmane.org/gmane.__science.biology.informatics.__conductor > <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 <mailto:hpages at="" fhcrc.org=""> > Phone: (206) 667-5791 <tel:%28206%29%20667-5791> > Fax: (206) 667-1319 <tel:%28206%29%20667-1319> > > -- 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
ADD REPLY •
0
Entering edit mode
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
ADD REPLY •

Login before adding your answer.

Similar Posts
Loading Similar Posts
Traffic: 725 users visited in the last hour
Content Search
Users
Tags
Badges
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6