Using GenomicRanges with table data
1
0
Entering edit mode
Tom Oates ▴ 60
@tom-oates-5703
Last seen 6.7 years ago
Hi I am trying to use GenomicRanges as part of an anlalysis of sequencing data. I have a number of files which I wish to use to make GRanges objects. For example: chr1 579578 579804 CpG_12 chr1 630418 630623 CpG_11 chr1 804552 804763 CpG_9 chr1 1307051 1307362 CpG_16 chr1 1323599 1323808 CpG_9 chr1 1350549 1350758 CpG_12 chr1 1403287 1403637 CpG_20 chr1 1418906 1419488 CpG_28 This file is sorted such that chr1 is followed chr2, 3, 4, 5 etc to 20 (as opposed to chr10, 11...19, 2, 3 etc) I use to make the GRanges object cpgi_gr<-GRanges(seqnames=Rle(cpgi$V1), ranges=IRanges(start=cpgi$V2,end=cpgi$V3), UCSC_AL_ID=cpgi$V4) but then if I examine seqnames(cpgi_gr) I get factor-Rle of length 89611 with 21 runs Lengths: 8952 5602 6133 4973 5840 4260 ... 3132 3607 2793 3175 3842 1419 Values : chr1 chr2 chr3 chr4 chr5 chr6 ... chr16 chr17 chr18 chr19 chr20 chrX Levels(21): chr1 chr10 chr11 chr12 chr13 chr14 chr15 ... chr5 chr6 chr7 chr8 chr9 chrX So the Values & Levels are not matching. I hope to give the GRanges object seqlengths of the chr lengths in the genome so I can then perform flank etc tasks on the data so it is crucial that the values & lengths match. I imagine that this problem is based around my not understanding either IRanges or Rle sufficiently but I have read help on Rle objects & IRanges and can't work out how to ensure that the formation of the GRanges object leads to the chr values matching Thanks
IRanges GenomicRanges IRanges GenomicRanges • 1.6k views
ADD COMMENT
0
Entering edit mode
@steve-lianoglou-2771
Last seen 13 months ago
United States
Hi On Friday, January 11, 2013, Tom Oates wrote: > Hi > I am trying to use GenomicRanges as part of an anlalysis of sequencing > data. > I have a number of files which I wish to use to make GRanges objects. > For example: > > chr1 579578 579804 CpG_12 > chr1 630418 630623 CpG_11 > chr1 804552 804763 CpG_9 > chr1 1307051 1307362 CpG_16 > chr1 1323599 1323808 CpG_9 > chr1 1350549 1350758 CpG_12 > chr1 1403287 1403637 CpG_20 > chr1 1418906 1419488 CpG_28 > > This file is sorted such that chr1 is followed chr2, 3, 4, 5 etc to 20 > (as opposed to chr10, 11...19, 2, 3 etc) > > I use to make the GRanges object > cpgi_gr<-GRanges(seqnames=Rle(cpgi$V1), > ranges=IRanges(start=cpgi$V2,end=cpgi$V3), > UCSC_AL_ID=cpgi$V4) > > but then if I examine > > seqnames(cpgi_gr) > > > I get > factor-Rle of length 89611 with 21 runs > Lengths: 8952 5602 6133 4973 5840 4260 ... 3132 3607 2793 > 3175 3842 1419 > Values : chr1 chr2 chr3 chr4 chr5 chr6 ... chr16 chr17 chr18 > chr19 chr20 chrX > Levels(21): chr1 chr10 chr11 chr12 chr13 chr14 chr15 ... chr5 chr6 > chr7 chr8 chr9 chrX > > > So the Values & Levels are not matching. I'm not sure what you mean by this, but your object looks more or less fine to me. If by the "levels not matching" you mean that chr10 comes (instead of chr2) after chr1 in a call to "sea levels", then that should really affect your analysis ... HTH, Steve I hope to give the GRanges > object seqlengths of the chr lengths in the genome so I can then > perform flank etc tasks on the data so it is crucial that the values & > lengths match. I imagine that this problem is based around my not > understanding either IRanges or Rle sufficiently but I have read help > on Rle objects & IRanges and can't work out how to ensure that the > formation of the GRanges object leads to the chr values matching > > Thanks > > _______________________________________________ > Bioconductor mailing list > Bioconductor@r-project.org <javascript:;> > https://stat.ethz.ch/mailman/listinfo/bioconductor > Search the archives: > http://news.gmane.org/gmane.science.biology.informatics.conductor > -- 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 [[alternative HTML version deleted]]
ADD COMMENT
0
Entering edit mode
Sorry for not being clear That is indeed what I mean by "levels not matching" So is it true that the chr1..chr2 etc as it appears in Levels will have no effect on further analysis? even though Values & Levels for seqnames do not match? On Fri, Jan 11, 2013 at 1:54 PM, Steve Lianoglou <mailinglist.honeypot at="" gmail.com=""> wrote: > Hi > > On Friday, January 11, 2013, Tom Oates wrote: >> >> Hi >> I am trying to use GenomicRanges as part of an anlalysis of sequencing >> data. >> I have a number of files which I wish to use to make GRanges objects. >> For example: >> >> chr1 579578 579804 CpG_12 >> chr1 630418 630623 CpG_11 >> chr1 804552 804763 CpG_9 >> chr1 1307051 1307362 CpG_16 >> chr1 1323599 1323808 CpG_9 >> chr1 1350549 1350758 CpG_12 >> chr1 1403287 1403637 CpG_20 >> chr1 1418906 1419488 CpG_28 >> >> This file is sorted such that chr1 is followed chr2, 3, 4, 5 etc to 20 >> (as opposed to chr10, 11...19, 2, 3 etc) >> >> I use to make the GRanges object >> cpgi_gr<-GRanges(seqnames=Rle(cpgi$V1), >> ranges=IRanges(start=cpgi$V2,end=cpgi$V3), >> UCSC_AL_ID=cpgi$V4) >> >> but then if I examine >> >> seqnames(cpgi_gr) >> >> >> I get >> factor-Rle of length 89611 with 21 runs >> Lengths: 8952 5602 6133 4973 5840 4260 ... 3132 3607 2793 >> 3175 3842 1419 >> Values : chr1 chr2 chr3 chr4 chr5 chr6 ... chr16 chr17 chr18 >> chr19 chr20 chrX >> Levels(21): chr1 chr10 chr11 chr12 chr13 chr14 chr15 ... chr5 chr6 >> chr7 chr8 chr9 chrX >> >> >> So the Values & Levels are not matching. > > > I'm not sure what you mean by this, but your object looks more or less fine > to me. > > If by the "levels not matching" you mean that chr10 comes (instead of chr2) > after chr1 in a call to "sea levels", then that should really affect your > analysis ... > > > HTH, > Steve > >> I hope to give the GRanges >> object seqlengths of the chr lengths in the genome so I can then >> perform flank etc tasks on the data so it is crucial that the values & >> lengths match. I imagine that this problem is based around my not >> understanding either IRanges or Rle sufficiently but I have read help >> on Rle objects & IRanges and can't work out how to ensure that the >> formation of the GRanges object leads to the chr values matching >> >> Thanks >> >> _______________________________________________ >> 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 > > > > -- > 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
Ideally your seqlevels would be in their natural ordering, i.e., chr1, chr2, etc. You could for example get a nicely ordered Seqinfo object from a BSgenome package or similar. For example seqinfo(gr) <- seqinfo(Hsapiens) But it doesn't matter as far as the integrity of your data. Michael On Fri, Jan 11, 2013 at 6:37 AM, Tom Oates <toates19@gmail.com> wrote: > Sorry for not being clear > That is indeed what I mean by "levels not matching" > So is it true that the chr1..chr2 etc as it appears in Levels will > have no effect on further analysis? even though Values & Levels for > seqnames do not match? > > On Fri, Jan 11, 2013 at 1:54 PM, Steve Lianoglou > <mailinglist.honeypot@gmail.com> wrote: > > Hi > > > > On Friday, January 11, 2013, Tom Oates wrote: > >> > >> Hi > >> I am trying to use GenomicRanges as part of an anlalysis of sequencing > >> data. > >> I have a number of files which I wish to use to make GRanges objects. > >> For example: > >> > >> chr1 579578 579804 CpG_12 > >> chr1 630418 630623 CpG_11 > >> chr1 804552 804763 CpG_9 > >> chr1 1307051 1307362 CpG_16 > >> chr1 1323599 1323808 CpG_9 > >> chr1 1350549 1350758 CpG_12 > >> chr1 1403287 1403637 CpG_20 > >> chr1 1418906 1419488 CpG_28 > >> > >> This file is sorted such that chr1 is followed chr2, 3, 4, 5 etc to 20 > >> (as opposed to chr10, 11...19, 2, 3 etc) > >> > >> I use to make the GRanges object > >> cpgi_gr<-GRanges(seqnames=Rle(cpgi$V1), > >> ranges=IRanges(start=cpgi$V2,end=cpgi$V3), > >> UCSC_AL_ID=cpgi$V4) > >> > >> but then if I examine > >> > >> seqnames(cpgi_gr) > >> > >> > >> I get > >> factor-Rle of length 89611 with 21 runs > >> Lengths: 8952 5602 6133 4973 5840 4260 ... 3132 3607 2793 > >> 3175 3842 1419 > >> Values : chr1 chr2 chr3 chr4 chr5 chr6 ... chr16 chr17 chr18 > >> chr19 chr20 chrX > >> Levels(21): chr1 chr10 chr11 chr12 chr13 chr14 chr15 ... chr5 chr6 > >> chr7 chr8 chr9 chrX > >> > >> > >> So the Values & Levels are not matching. > > > > > > I'm not sure what you mean by this, but your object looks more or less > fine > > to me. > > > > If by the "levels not matching" you mean that chr10 comes (instead of > chr2) > > after chr1 in a call to "sea levels", then that should really affect your > > analysis ... > > > > > > HTH, > > Steve > > > >> I hope to give the GRanges > >> object seqlengths of the chr lengths in the genome so I can then > >> perform flank etc tasks on the data so it is crucial that the values & > >> lengths match. I imagine that this problem is based around my not > >> understanding either IRanges or Rle sufficiently but I have read help > >> on Rle objects & IRanges and can't work out how to ensure that the > >> formation of the GRanges object leads to the chr values matching > >> > >> Thanks > >> > >> _______________________________________________ > >> 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 > > > > > > > > -- > > 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
one day I found myself in this situation and wrote the Stupid Man's version of the above... i.e. "just do what I mean" (yes that's bad) addSeqinfo <- function (x) { require(rtracklayer) ## and assumes that genome(x) has been set stopifnot(is(x, "SummarizedExperiment") || is(x, "GenomicRanges")) gen <- ifelse(is(x, "SummarizedExperiment"), genome(rowData(x)), genome(x)) gen <- unique(gen) if is.na(gen)) stop("You must assign genome(x) first") if (is(x, "SummarizedExperiment")) { seqinfo(rowData(x)) <- SeqinfoForBSGenome(gen)[seqlevels(rowData(x))] } else { seqinfo(x) <- SeqinfoForBSGenome(gen)[seqlevels(x)] } return(x) } Granted it is not beautiful but when I automatically pull some HuEx, methylation array, eRRBS or summarized RNAseq data from GEO it can be nice for the whole chain to just "know what to do" based on the genome assembly of the annotations used to create the GRanges or SE: R> unique(genome(BcellALL)) ## [1] "hg19" R> head(as.data.frame(seqinfo(BcellALL))) ## data frame with 6 rows and 3 columns ## seqlengths isCircular genome ## <integer> <logical> <character> ## chr1 NA NA hg19 ## ... R> BcellALL <- addSeqinfo(BcellALL) R> head(as.data.frame(seqinfo(BcellALL))) ## data frame with 6 rows and 3 columns ## seqlengths isCircular genome ## <integer> <logical> <character> ## chr1 249250621 FALSE hg19 ## ... Needless to say I am a big fan of BSgenome and rtracklayer -- thank you Herve and Michael for putting them together (and yes I will get around to testing import.bigBed() when I can... I cheated and used bigBedToBed for the task originally motivating the patch, however). On Fri, Jan 11, 2013 at 7:14 AM, Michael Lawrence <lawrence.michael@gene.com> wrote: > Ideally your seqlevels would be in their natural ordering, i.e., chr1, > chr2, etc. You could for example get a nicely ordered Seqinfo object from a > BSgenome package or similar. > > For example > seqinfo(gr) <- seqinfo(Hsapiens) > > But it doesn't matter as far as the integrity of your data. > > Michael > > > On Fri, Jan 11, 2013 at 6:37 AM, Tom Oates <toates19@gmail.com> wrote: > > > Sorry for not being clear > > That is indeed what I mean by "levels not matching" > > So is it true that the chr1..chr2 etc as it appears in Levels will > > have no effect on further analysis? even though Values & Levels for > > seqnames do not match? > > > > On Fri, Jan 11, 2013 at 1:54 PM, Steve Lianoglou > > <mailinglist.honeypot@gmail.com> wrote: > > > Hi > > > > > > On Friday, January 11, 2013, Tom Oates wrote: > > >> > > >> Hi > > >> I am trying to use GenomicRanges as part of an anlalysis of sequencing > > >> data. > > >> I have a number of files which I wish to use to make GRanges objects. > > >> For example: > > >> > > >> chr1 579578 579804 CpG_12 > > >> chr1 630418 630623 CpG_11 > > >> chr1 804552 804763 CpG_9 > > >> chr1 1307051 1307362 CpG_16 > > >> chr1 1323599 1323808 CpG_9 > > >> chr1 1350549 1350758 CpG_12 > > >> chr1 1403287 1403637 CpG_20 > > >> chr1 1418906 1419488 CpG_28 > > >> > > >> This file is sorted such that chr1 is followed chr2, 3, 4, 5 etc to 20 > > >> (as opposed to chr10, 11...19, 2, 3 etc) > > >> > > >> I use to make the GRanges object > > >> cpgi_gr<-GRanges(seqnames=Rle(cpgi$V1), > > >> ranges=IRanges(start=cpgi$V2,end=cpgi$V3), > > >> UCSC_AL_ID=cpgi$V4) > > >> > > >> but then if I examine > > >> > > >> seqnames(cpgi_gr) > > >> > > >> > > >> I get > > >> factor-Rle of length 89611 with 21 runs > > >> Lengths: 8952 5602 6133 4973 5840 4260 ... 3132 3607 2793 > > >> 3175 3842 1419 > > >> Values : chr1 chr2 chr3 chr4 chr5 chr6 ... chr16 chr17 chr18 > > >> chr19 chr20 chrX > > >> Levels(21): chr1 chr10 chr11 chr12 chr13 chr14 chr15 ... chr5 chr6 > > >> chr7 chr8 chr9 chrX > > >> > > >> > > >> So the Values & Levels are not matching. > > > > > > > > > I'm not sure what you mean by this, but your object looks more or less > > fine > > > to me. > > > > > > If by the "levels not matching" you mean that chr10 comes (instead of > > chr2) > > > after chr1 in a call to "sea levels", then that should really affect > your > > > analysis ... > > > > > > > > > HTH, > > > Steve > > > > > >> I hope to give the GRanges > > >> object seqlengths of the chr lengths in the genome so I can then > > >> perform flank etc tasks on the data so it is crucial that the values & > > >> lengths match. I imagine that this problem is based around my not > > >> understanding either IRanges or Rle sufficiently but I have read help > > >> on Rle objects & IRanges and can't work out how to ensure that the > > >> formation of the GRanges object leads to the chr values matching > > >> > > >> Thanks > > >> > > >> _______________________________________________ > > >> 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 > > > > > > > > > > > > -- > > > 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]] > > _______________________________________________ > 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 > -- *A model is a lie that helps you see the truth.* * * Howard Skipper<http: cancerres.aacrjournals.org="" content="" 31="" 9="" 1173.full.pdf=""> [[alternative HTML version deleted]]
ADD REPLY

Login before adding your answer.

Traffic: 940 users visited in the last hour
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