GRs are a great data structure. But, standard bioinformatic file formats
(BED, WIG, BAM) don't always fit 1:1 with the "organic" beginnings of some
projects. The GenomicRanges infrastructure isn't on the radar of every R
developer, and some useful data can be found in ugly formats. Wouldn't it
be handy if users could easily turn those blobs of data into something that
export() can handle?
To better understand Michael's point of view, and as someone who has seen
firsthand the nontrivial amount of work required to maintain rtracklayer as
a high-performance import library, I wrote a few trivial Perl scripts to
convert some mouse data from assorted wacky tabular formats to standard
BED6 files. Besides, once I have BEDs, I can Tabix them, which speeds up
operations.
I noticed that, when I imported the resulting BED files, where I had cloned
the base position into identical "chromStart" and "chromEnd" coordinates,
the import.bed() function assumed that I meant for them to be UCSC-style,
and therefore gave everything negative widths. (On the bright side, this
also explained why liftOver wasn't doing anything useful with the results)
packageVersion('rtracklayer')
## [1] '1.21.12'
wacky <- import('converted.theirSillyFormat.mm9.bed.gz', genome='mm9')
mm9toHg19 <- import('mm9ToHg19.over.chain')
empty <- liftOver(wacky, mm9toHg19)
length(unlist(empty))/length(wacky)
## [1] 0
## cursing ensues
That's not so good. If I wasn't already aware of the insanity of
UCSC-style indexing, this could have been a problem in and of itself. (As
it was, I fixed it)
slow <- read.table('theirSillyFormat.mm9.txt.gz')
## time passes...
aGR <- df2GR(slow)
mm9toHg19 <- import('mm9ToHg19.over.chain')
full <- liftOver(slow, mm9toHg19)
length(unlist(full))/length(slow)
## [1] 0.549
## better late than never
So: turning nonstandard data into standard data with "standard" (grr) UCSC
assumptions took longer than simply brute-forcing the issue with
read.table().
After importing the files with the "appropriate" (chrom, base-1, base)
indexing, I then went to liftOver() a related GRanges. (Note that the
related GRanges was from a BED file submitted by guys at the Broad, so any
wacky formatting wasn't an issue in this case; here I wanted to control for
any other possible fubars).
foo <- import('an.RRBS.file.mm9.bed.gz', genome='mm9')
mm9toHg19 <- import('mm9ToHg19.over.chain')
bar <- liftOver(foo, mm9toHg19)
length(unlist(bar)) / length(foo)
## [1] 0.622
Needless to say, this was a hell of a lot faster than importing the
corresponding file as a table. However, for the wacky file formats, it was
more time & trouble to decode all the assumptions prior to liftOver() than
it would have been to use granges(read.table('wacky.file.format.csv.gz')).
Ugly and sad, but still true!
So, in conclusion, sometimes it might just be better to import one of those
wacky file formats to a data.frame D or what have you, and use granges(D).
Or, since I'm one of the people that wrote a df2GR() function, I just use
that.
--t
On Mon, Oct 7, 2013 at 9:23 AM, Steve Lianoglou <lianoglou.st...@gene.com>wrote:
> Hi,
>
> +1 from me, too ... I've also had a similar conversion function
> (data.[frame|table] <--> GRanges) in my toolbelt which I found quite
> useful over the years.
>
> -steve
>
> On Sun, Oct 6, 2013 at 5:00 PM, Kasper Daniel Hansen
> <kasperdanielhan...@gmail.com> wrote:
> > This is a convenience function, which provably has saved tons of time for
> > me and others. I get lots of data from various excel/cvs files lying
> > around various places, and these files _always_ have a clear path to a
> > GRanges. Perhaps you never have to deal with this kind of data, but we
> are
> > a few experienced users who find it extremely handy and would like it to
> be
> > in a more centralized place.
> >
> >
> >
> >
> >
> > On Sun, Oct 6, 2013 at 4:26 PM, Michael Lawrence
> > <lawrence.mich...@gene.com>wrote:
> >
> >> I'm still unconvinced that there is an obvious, general path from
> >> data.frame -> GRanges. It's usually easy enough to just call GRanges(),
> >> often of the pattern with(df, GRanges(...)). Moreover, it's unusual for
> me
> >> to encounter genomic data in data.frames.
> >>
> >>
> >>
> >>
> >> On Sun, Oct 6, 2013 at 8:37 AM, Kasper Daniel Hansen <
> >> kasperdanielhan...@gmail.com> wrote:
> >>
> >>> Also, it goes without saying that I am happy to provide a patch for
> >>> GenomicRanges, and check that it passes R CMD check, to minimize the
> work
> >>> of the maintainer.
> >>>
> >>> Kasper
> >>>
> >>>
> >>> On Sun, Oct 6, 2013 at 9:28 AM, Kasper Daniel Hansen <
> >>> kasperdanielhan...@gmail.com> wrote:
> >>>
> >>> > bsseq::data.frame2GRanges does the obvious step of converting a
> >>> data.frame
> >>> > to GRanges. It has a couple of bells and whistles where strand can
> be
> >>> > ignored and additional columns (apart from genomic location) may be
> >>> ignore
> >>> > in the output object.
> >>> >
> >>> > I (and now quite a few other people) use this function almost every
> day.
> >>> > I have seen other implementations in other packages, suggesting
> this is
> >>> > not just something I (we) use.
> >>> >
> >>> > I suggests adding this function to GenomicRanges. I am happy to
> support
> >>> > it going forward.
> >>> >
> >>> > Using this function we could also add an as(x, "GRanges") method for
> >>> > x=data.frame, but I still suggest keeping the basic function for the
> >>> > extended functionality it provides.
> >>> >
> >>> > Best,
> >>> > Kasper
> >>> >
> >>>
> >>> [[alternative HTML version deleted]]
> >>>
> >>> _______________________________________________
> >>> Bioc-devel@r-project.org mailing list
> >>> https://stat.ethz.ch/mailman/listinfo/bioc-devel
> >>>
> >>
> >>
> >
> > [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > Bioc-devel@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/bioc-devel
>
>
>
> --
> Steve Lianoglou
> Computational Biologist
> Bioinformatics and Computational Biology
> Genentech
>
> _______________________________________________
> Bioc-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/bioc-devel
>
--
*He that would live in peace and at ease, *
*Must not speak all he knows, nor judge all he sees.*
*
*
Benjamin Franklin, Poor Richard's
Almanack<http://archive.org/details/poorrichardsalma00franrich>
[[alternative HTML version deleted]]
_______________________________________________
Bioc-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/bioc-devel
