[BioC] some random comments on GRanges

[Kasper Daniel Hansen]

What follows is some random comments on the use of GRanges for storing
data in the elementMetadata slot.  My use case is a lot of very small
ranges (>10M ranges of width 1), which may be a special border case
(but very important to me :).  I am particular interested in fast
subsetting like subsetByOverlaps.

Comments:

* I find column subsetting on the elementMetadata to be really slow,
slower than it need to be I would say (I suspect some validity
checking, see comment below on row subsetting).

* There is no easy shortcut to get a specific metadata column (gr[,
"score"] returns a GRange).  I would find it natural that gr$score
would do that.  This is important in order to get at the data stores
in the GRange

* while there is a as.data.frame, as(from, "data.frame") does not work
and I would also appreciate an as(from, "GRanges") method with from
being a data.frame.  I have a quick function at the bottom of my
email, for doing this.

* Since start, end, width are reserved column names for GRanges, it
would be nice if the GRanges constructor was extended to allow
  GRanges(seqnames = "chr1", start = 1, end = 10)
instead of now where we need
  GRanges(seqnames = "chr1", ranges = IRanges(start = 1, end = 10))
(ok, this is pretty minor, but it is convenient)

* There is no sort method implemented for GRanges; it returns an error

* is.unsorted(xx) seems to always return FALSE with a warning when xx
is an IRange or a GRange

*It seems to me that row subsetting, specifically subsetByOverlaps is
a bit too slow.  The find overlaps part is blazingly fast (when ranges
is a single IRange); it seems like the initialization and subsetting
part takes "more time than necessary".  This is based on comparing a
(long) GRanges with around 8 elementMetadata columns to another
approach where instead of using elementMetadata I use two elements of
a list: a GRanges without elementMetadata and where the
elementMetadata is stored in a separate data.frame - I then use
findOverlaps to get me the indices (when I compare using a GRanges to
using a GRanges and separate matrix, I use the "position" GRanges to
get my indexes, which I use to subset the matrix with).

Specifically here gr is a GRanges with length 25M (no elementMetadata)
and mat is a 25M x 8 matrix

grIn <- gr
elementMetadata(grIn) <- mat

## The "classic" approach:
system.time({
grOut1 <- subsetByOverlaps(grIn, GRanges(seqnames = "chr1", ranges =
IRanges(start = 1, end = 10^7)))
})  # 4.7 secs

## Keeping gr and mat separate:
system.time({
    qh <- queryHits(findOverlaps(grIn, GRanges(seqnames = "chr1",
ranges = IRanges(start = 1, end = 10^7))))
    grOut2 <- list(gr[qh],  mat = mat[qh,])
}) #  0.9 secs

## Doing the above, but forming a GRange after the subsetting adds
just a small overhead
system.time({
    qh <- queryHits(findOverlaps(grIn, GRanges(seqnames = "chr1",
ranges = IRanges(start = 1, end = 10^7))))
    grOut3 = gr[qh]
    elementMetadata(grOut3) <- mat[qh,]
}) # 1.1 secs


* For really large query GRanges timing results (not shown) seems to
indicate that the line
    query[!is.na(findOverlaps(query, subject, maxgap = maxgap,
        minoverlap = minoverlap, type = type, select = "first"))]
in subsetByOverlaps(query = "GRanges", subject = "GRanges") could be
made faster by changing it to
    query[queryHits(findOverlaps(query, subject, maxgap = maxgap,
        minoverlap = minoverlap, type = type))]
(my timing is for about a query length of 30M and a subject length of 1)

* I find it counter-intuitive that the following code doesn't work:
R> example(GRanges)
R> seqnames(gr) = "chr1"
or perhaps more intuitive
R> gr1 = gr[seqnames(gr) == "Chrom1"]
R> seqnames(gr1) = "chr1"
In general I don't like that seqnames is a factor.  This is fine when
I just consider a single GRange, but when I put together multiple
GRanges from different sources, I sometimes get into problems.
Specifically, in human, it is not always true that the various contigs
(like "chr10_random") are part of the data sources.  I would have
preferred it to be a character and then have seqlengths essentially be
NA in case a character does not appear in the names(seqlengths).  I
know that I have complained about this a long time ago and it is
probably too late to change, but still :)


Kasper

data.frame2GRanges <- function(object, keepColumns = TRUE) {
    stopifnot(class(object) == "data.frame")
    stopifnot(all(c("chr", "start", "end") %in% names(object)))
    if("strand" %in% names(object)) {
        if(is.numeric(object$strand)) {
            strand <- ifelse(object$strand == 1, "+", "*")
            strand[object$strand == -1] <- "-"
            object$strand <- strand
        }
        gr <- GRanges(seqnames = object$chr,
                      ranges = IRanges(start = object$start, end = object$end),
                      strand = object$strand)
    } else {
        gr <- GRanges(seqnames = object$chr,
                      ranges = IRanges(start = object$start, end = object$end))
    }
    if(keepColumns) {
        dt <- as(object[, setdiff(names(object), c("chr", "start",
"end", "strand"))],
                     "DataFrame")
        elementMetadata(gr) <- dt
    }
    names(gr) <- rownames(object)
    gr
}