R/AnnotationTrack-class.R

In ivanek/Gviz: Plotting data and annotation information along genomic coordinates

#' @include DataTrack-class.R
#' @include ReferenceTrack-class.R

NULL

## AnnotationTrack Class -----------------------------------------------------

#' AnnotationTrack class and methods
#'
#' A fairly generic track object for arbitrary genomic range annotations, with
#' the option of grouped track items. The extended `DetailsAnnotationTrack`
#' provides a more flexible interface to add user-defined custom
#' information for each range.
#'
#' @name AnnotationTrack-class
#'
#' @template GdObject-class_slot
#' @template RangeTrack-class_slot
#' @template StackedTrack_slot
#' @template AnnotationTrack_slot
#'
#' @template AnnotationTrack_param
#'
#' @author Florian Hahne, Arne Mueller
#'
#' @inherit GdObject-class seealso
#'
#' @examples
#' ## An empty object
#' AnnotationTrack()
#'
#' ## Construct from individual arguments
#' st <- c(2000000, 2070000, 2100000, 2160000)
#' ed <- c(2050000, 2130000, 2150000, 2170000)
#' str <- c("-", "+", "-", "-")
#' gr <- c("Group1", "Group2", "Group1", "Group3")
#'
#' annTrack <- AnnotationTrack(
#'     start = st, end = ed, strand = str, chromosome = 7,
#'     genome = "hg19", feature = "test", group = gr,
#'     id = paste("annTrack item", 1:4),
#'     name = "generic annotation", stacking = "squish"
#' )
#'
#' ## Or from a data.frame
#' df <- data.frame(
#'     start = st, end = ed, strand = str, id = paste("annTrack item", 1:4),
#'     feature = "test", group = gr
#' )
#' annTrack <- AnnotationTrack(
#'     range = df, genome = "hg19", chromosome = 7,
#'     name = "generic annotation", stacking = "squish"
#' )
#'
#' ## Or from a GRanges object
#' gr <- GRanges(
#'     seqnames = "chr7", range = IRanges(start = df$start, end = df$end),
#'     strand = str
#' )
#' genome(gr) <- "hg19"
#' mcols(gr) <- df[, -(1:3)]
#' annTrack <- AnnotationTrack(
#'     range = gr, name = "generic annotation",
#'     stacking = "squish"
#' )
#'
#' ## Finally from a GRangesList
#' grl <- split(gr, values(gr)$group)
#' AnnotationTrack(grl)
#' \dontshow{
#' ## For some annoying reason the postscript device does not know about
#' ## the sans font
#' if (!interactive()) {
#'     font <- ps.options()$family
#'     displayPars(annTrack) <- list(fontfamily = font, fontfamily.title = font)
#' }
#' }
#'
#' ## Plotting
#' plotTracks(annTrack)
#'
#' ## Track names
#' names(annTrack)
#' names(annTrack) <- "foo"
#' plotTracks(annTrack)
#'
#' ## Subsetting and splitting
#' subTrack <- subset(annTrack, to = 2155000)
#' length(subTrack)
#' subTrack[1:2]
#' split(annTrack, c(1, 2, 1, 2))
#'
#' ## Accessors
#' start(annTrack)
#' end(annTrack)
#' width(annTrack)
#' position(annTrack)
#' width(subTrack) <- width(subTrack) + 1000
#'
#' strand(annTrack)
#' strand(subTrack) <- "-"
#'
#' chromosome(annTrack)
#' chromosome(subTrack) <- "chrX"
#'
#' genome(annTrack)
#' genome(subTrack) <- "mm9"
#'
#' range(annTrack)
#' ranges(annTrack)
#'
#' ## Annotation
#' identifier(annTrack)
#' identifier(annTrack, "lowest")
#' identifier(subTrack) <- "bar"
#'
#' feature(annTrack)
#' feature(subTrack) <- "foo"
#'
#' values(annTrack)
#'
#' ## Grouping
#' group(annTrack)
#' group(subTrack) <- "Group 1"
#' chromosome(subTrack) <- "chr7"
#' plotTracks(subTrack)
#'
#' ## Stacking
#' stacking(annTrack)
#' stacking(annTrack) <- "dense"
#' plotTracks(annTrack)
#'
#' ## coercion
#' as(annTrack, "data.frame")
#' as(annTrack, "UCSCData")
#'
#' ## HTML image map
#' coords(annTrack)
#' tags(annTrack)
#' annTrack <- plotTracks(annTrack)$foo
#' coords(annTrack)
#' tags(annTrack)
#'
#' ## DetailsAnnotationTrack
#' library(lattice) # need to use grid grapics
#'
#' ## generate two random distributions per row (probe/feature)
#' ## the difference between the distributions increases from probe 1 to 4
#' m <- matrix(c(rgamma(400, 1)), ncol = 100)
#' m[, 51:100] <- m[, 51:100] + 0:3
#' ## rownames must be accessible by AnnotationTrack element identifier
#' rownames(m) <- identifier(annTrack, "lowest")
#'
#' ## create a lattice density plot for the values (signals) of the two groups
#' ## as the chart must be placed into a pre-set grid view port we have to use
#' ## print without calling plot.new! Note, use a common prefix for all lattice.
#' ## Avoid wasting space by removing y-axis decorations.
#'
#' ## Note, in this example 'm' will be found in the environment the 'details'
#' ## function is defined in. To avoid overwriting 'm' you should use a closure
#' ## or environment to access 'm'.
#' details <- function(identifier, ...) {
#'     d <- data.frame(signal = m[identifier, ], group = rep(c("grp1", "grp2"), each = 50))
#'     print(densityplot(~signal,
#'         group = group, data = d, main = identifier,
#'         scales = list(draw = FALSE, x = list(draw = TRUE)), ylab = "", xlab = "",
#'     ), newpage = FALSE, prefix = "plot")
#' }
#'
#' deTrack <- AnnotationTrack(
#'     range = gr, genome = "hg19", chromosome = 7,
#'     name = "generic annotation with details per entry", stacking = "squish",
#'     fun = details, details.ratio = 1
#' )
#'
#' plotTracks(deTrack)
#'
#' set.seed(1234)
#' deTrack <- AnnotationTrack(
#'     range = gr, genome = "hg19", chromosome = 7,
#'     name = "generic annotation with details per entry",
#'     stacking = "squish", fun = details,
#'     details.ratio = 1, selectFun = function(...) {
#'         sample(c(FALSE, TRUE), 1)
#'     }
#' )
#'
#' plotTracks(deTrack)
#' @importClassesFrom rtracklayer UCSCData
#'
#' @exportClass AnnotationTrack
setClass("AnnotationTrack",
    contains = "StackedTrack",
    prototype = prototype(
        columns = c("feature", "group", "id"),
        stacking = "squish",
        name = "AnnotationTrack",
        dp = DisplayPars(
            arrowHeadWidth = 30,
            arrowHeadMaxWidth = 40,
            cex.group = 0.6,
            cex = 1,
            col.line = "darkgray",
            col = "transparent",
            featureAnnotation = NULL,
            fill = "lightblue",
            fontfamily.group = "sans",
            fontcolor.group = .DEFAULT_SHADED_COL,
            fontcolor.item = "white",
            fontface.group = 2,
            fontsize.group = 12,
            groupAnnotation = NULL,
            just.group = "left",
            lex = 1,
            lineheight = 1,
            lty = "solid",
            lwd = 1,
            mergeGroups = FALSE,
            min.height = 3,
            min.width = 1,
            rotation = 0,
            rotation.group = 0,
            rotation.item = 0,
            shape = "arrow",
            showFeatureId = FALSE,
            showId = FALSE,
            showOverplotting = FALSE,
            size = 1
        )
    )
)

## Initialize ----------------------------------------------------------------

## Essentially we just check for the correct GRanges columns here
#' @describeIn AnnotationTrack-class Show method.
#' @export
setMethod("initialize", "AnnotationTrack", function(.Object, ...) {
    if (is.null(list(...)$range) && is.null(list(...)$genome) && is.null(list(...)$chromosome)) {
        return(.Object)
    }
    ## the display parameter defaults
    .makeParMapping()
    .Object <- .updatePars(.Object, "AnnotationTrack")
    range <- list(...)$range
    if (!is.null(range) && length(.Object)) {
        if (!all(.Object@columns %in% colnames(values(range)))) {
            stop(
                "Problem initializing AnnotationTrack, need the following columns:",
                toString(.Object@columns, collapse = ", ")
            )
        }
        grp <- if (is(.Object, "GeneRegionTrack")) values(range)$transcript else values(range)$group
        if (any(vapply(split(as.character(strand(range)), grp), function(x) length(unique(x)), numeric(1)) != 1)) {
            stop("Grouped elments of a RangeTrack can not be on opposing strands")
        }
    }
    .Object <- callNextMethod()
    return(.Object)
})

## The file-based version of the AnnotationTrack class. This will mainly provide a means to dispatch to
## a special 'subset' method which should stream the necessary data from disk.

#' @describeIn AnnotationTrack-class The file-based version of the `AnnotationTrack-class`.
#' @exportClass ReferenceAnnotationTrack
setClass("ReferenceAnnotationTrack", contains = c("AnnotationTrack", "ReferenceTrack"))

## Initialize ----------------------------------------------------------------

## This just needs to set the appropriate slots that are being inherited from ReferenceTrack because the
## multiple inheritance has some strange features with regards to method selection
#' @describeIn AnnotationTrack-class Initialize.
#' @export
setMethod("initialize", "ReferenceAnnotationTrack", function(.Object, stream, reference, mapping = list(),
                                                             args = list(), defaults = list(), ...) {
    .Object <- selectMethod("initialize", "ReferenceTrack")(.Object = .Object, reference = reference, stream = stream,
        mapping = mapping, args = args, defaults = defaults)
    .Object <- callNextMethod(.Object, ...)
    return(.Object)
})

## Constructor ---------------------------------------------------------------

## Constructor. The following arguments are supported:
##    o range: a data.frame or a GRanges object containing the information
##       about the track items. If a data.frame, it needs to be coerceable
##     to a GRanges object, i.e., it needs at least the mandatory 'start', 'stop' and
##     'strand' columns. Additional optional columns are:
##        - feature: the type of the item. Can be mapped to colors via the DisplayPars.
##        - group: a grouping factor to connect track items.
##        - id: a unique identifier for a feature. This will be plotted if showId==TRUE
##             Note that internally we use the value of ID as the seqnames slot in the
##             internal GRanges object. Defaults for all missing columns are generated.
##    Instead of using the 'range' parameter, all these values can also be passed as
##    individual vectors, in which case they need to be of similar length.
##    o start, end, width: numeric vectors of the item start and end coordinates
##    o strand: the strand information may be provided in the form '+' for
##       the Watson strand, '-' for the Crick strand or '*' for any of the two.
##    o feature, group, id: individual vectors of equal length as described above.
##    o genome, chromosome: the reference genome and active chromosome for the track.
##    o stacking: character controlling the stacking of overlapping items. One in 'hide',
##       'dense', 'squish', 'pack' or 'full'.
##    o name: the name of the track. This will be used for the title panel.
## All additional items in ... are being treated as further DisplayParameters

#' @describeIn AnnotationTrack-class Constructor function for
#' `AnnotationTrack-class`

#' @return
#' The return value of the constructor function is a new object of class
#' \code{AnnotationTrack} or of class \code{DetailsAnnotationTrack}, depending
#' on the constructor arguments. Typically the user will not have to be
#' troubled with this distinction and can rely on the constructor to make the
#' right choice.
#'
#' @section Objects from the class:
#'
#' Objects can be created using the constructor function
#' \code{AnnotationTrack}.
#' @export
AnnotationTrack <- function(range = NULL, start = NULL, end = NULL, width = NULL, feature, group, id, strand, chromosome,
                            genome, stacking = "squish", name = "AnnotationTrack", fun, selectFun, importFunction,
                            stream = FALSE, ...) {
    ## Some defaults
    covars <- .getCovars(range)
    isStream <- FALSE
    if (!is.character(range)) {
        n <- max(c(length(start), length(end), length(width)), nrow(covars))
        if (is.null(covars[["feature"]]) && missing(feature)) {
            feature <- rep("unknown", n)
        }
        if (is.null(covars[["id"]]) && missing(id)) {
            id <- make.unique(rep(if (!is.null(feature)) as.character(feature) else covars[["feature"]], n)[seq_len(n)])
        }
        if (is.null(covars[["group"]]) && missing(group)) {
            group <- seq_len(n)
        }
    }
    ## Build a GRanges object from the inputs
    .missingToNull(c("feature", "group", "id", "strand", "chromosome", "importFunction", "genome"))
    args <- list(feature = feature, group = group, id = id, strand = strand, chromosome = chromosome, genome = genome)
    defs <- list(feature = "unknown", group = "unknown", id = "unknown", strand = "*", density = 1, chromosome = "chrNA", genome = NA)
    range <- .buildRange(
        range = range, groupId = "group", start = start, end = end, width = width,
        args = args, defaults = defs, chromosome = chromosome, trackType = "AnnotationTrack",
        importFun = importFunction, stream = stream
    )
    if (is.list(range)) {
        isStream <- TRUE
        slist <- range
        range <- GRanges()
    }
    ## Pipes have a special meaning for merged groups, so we can't have them in the initial group vector
    mcols(range)[["group"]] <- gsub("|", "", mcols(range)[["group"]], fixed = TRUE)
    ## If no chromosome was explicitely asked for we just take the first one in the GRanges object
    if (missing(chromosome) || is.null(chromosome)) {
        chromosome <- if (length(range) > 0) .chrName(as.character(seqnames(range)[1])) else "chrNA"
    }
    ## And finally the object instantiation, we have to distinguish between DetailsAnnotationTracks and normal ones
    genome <- .getGenomeFromGRange(range, ifelse(is.null(genome), character(), genome[1]))
    if (missing(fun)) {
        if (!isStream) {
            return(new("AnnotationTrack",
                chromosome = as.character(chromosome[1]), range = range,
                name = name, genome = genome, stacking = stacking, ...
            ))
        } else {
            ## A bit hackish but for some functions we may want to know which track type we need but at the
            ## same time we do not want to enforce this as an additional argument
            e <- new.env()
            e[["._trackType"]] <- "AnnotationTrack"
            environment(slist[["stream"]]) <- e
            return(new("ReferenceAnnotationTrack",
                chromosome = as.character(chromosome[1]), range = range,
                name = name, genome = genome, stacking = stacking, stream = slist[["stream"]], reference = slist[["reference"]],
                mapping = slist[["mapping"]], args = args, defaults = defs, ...
            ))
        }
    } else {
        if (!is.function(fun)) {
            stop("'fun' must be a function")
        }
        if (missing(selectFun)) {
            selectFun <- function(...) {
                return(TRUE)
            }
        }
        if (!is.function(selectFun)) {
            stop("'selectFun' must be a function")
        }
        return(new("DetailsAnnotationTrack",
            chromosome = as.character(chromosome[1]), range = range,
            name = name, genome = genome, stacking = stacking, fun = fun, selectFun = selectFun, ...
        ))
    }
}

#' @describeIn AnnotationTrack-class directly extends `AnnotationTrack.`
#' @exportClass DetailsAnnotationTrack
setClass("DetailsAnnotationTrack",
    contains = "AnnotationTrack",
    representation = representation(fun = "function", selectFun = "function"),
    prototype = prototype(
        fun = function(...) {},
        selectFun = function(...) {
            return(TRUE)
        },
        dp = DisplayPars(
            details.minWidth = 100,
            details.ratio = Inf,
            details.size = 0.5,
            detailsBorder.col = "darkgray",
            detailsBorder.fill = "transparent",
            detailsBorder.lty = "solid",
            detailsBorder.lwd = 1,
            detailsConnector.cex = 1,
            detailsConnector.col = "darkgray",
            detailsConnector.lty = "dashed",
            detailsConnector.lwd = 1,
            detailsConnector.pch = 20,
            detailsFunArgs = list(),
            groupDetails = FALSE
        )
    )
)

#' @describeIn AnnotationTrack-class Constructor function for
#' `DetailsAnnotationTrack-class`
#'
#' The `DetailsAnnotationTrack` class directly extends `AnnotationTrack.`
#' The purpose of this track type is to add an arbitrarily detailed plot
#' section (typically consisting of additional quantitative data) for each
#' range element of an `AnnotationTrack.` This allows a locus wide view of
#' annotation elements together with any kind of details per feature or element
#'  that may for instance provide insight on how some complex quantitative
#'  measurements change according to their position in a locus. If the
#'  quantitative data is too complex for a `DataTrack` e.g. because it requires
#'  extra space or a trellis-like representation, a `DetailsAnnotationTrack` can
#'  be used instead. Example: An `AnnotationTrack` shows the positions of a
#'  number of probes from a microarray, and you want a histogram of the signal
#'  intensity distribution derived from all samples at each of these probe
#'  location. Another example usage would be to show for each element of an
#'  `AnnotationTrack` an xy-plot of the signal against some clinical measurement
#'  such as blood pressure. The limitation for applications of this type of
#'  track is basically only the available space of the device you are
#'  plotting to.
#'
#' This flexibility is possible by utilizing a simple function model
#' to perform all the detailed plotting. The functionality of this plotting
#' function fun is totally up to the user, and the function environment is
#' prepared in a way that all necessary information about the plotted
#' annotation feature is available. To restrict the details section to only
#' selected number of annotation features one can supply another function
#' `selectFun`, which decides for each feature separately whether details are
#' available or not. Finally, an arbitrary number of additional arguments can
#' be passed on to these two function by means of the `detailsFunArgs` display
#' parameter. This is expected to be a named list, and all list elements are
#' passed along to the plotting function fun and to the selector function
#' `selectFun` as additional named arguments. Please note that some argument
#' names like `start`, `end` or `identifier` are reserved and can not be used
#' in the `detailsFunArgs` list. For examples of plotting functions,
#' see the 'Examples' section.
#' @export
DetailsAnnotationTrack <- function(...) AnnotationTrack(...)

## Initialize ----------------------------------------------------------------

#' @describeIn AnnotationTrack-class Initialize.
#' @export
setMethod("initialize", "DetailsAnnotationTrack", function(.Object, fun, selectFun, ...) {
    ## the diplay parameter defaults
    .Object <- .updatePars(.Object, "DetailsAnnotationTrack")
    .makeParMapping()
    .Object@fun <- fun
    .Object@selectFun <- selectFun
    .Object <- callNextMethod(.Object, ...)
    return(.Object)
})

## General accessors ---------------------------------------------------------
## Annotation Accessors ------------------------------------------------------

#' @describeIn AnnotationTrack-class extract the group membership for all track items.
#' @export
setMethod("group", "AnnotationTrack", function(GdObject) .getAnn(GdObject, "group"))

#' @describeIn AnnotationTrack-class replace the grouping information for track items.
#' The replacement value must be a factor of appropriate length or another
#' vector that can be coerced into such.
#' @export
setReplaceMethod("group", signature("AnnotationTrack", "character"), function(GdObject, value) .setAnn(GdObject, value, "group"))

## Context-dependent meta-accessors to the identifier data of an AnnotationTrack and a GeneRegionTrack. For the former, those will
## be the content of the 'id', the 'group' or the 'feature' metadata column, for the latter, one in
## the selection of 'symbol', 'gene', 'transcript', 'feature' or 'exon'. For historical reasons, logical values are also allowed, where TRUE
## is equivalent to 'symbol' and FALSE is equivalent to 'gene'. If not provided explicitely in the 'type' argument, the identifier
## type is inferred from the 'transcriptAnnotation' display parameter (a.k.a. 'geneSymbols') for a GeneRegionTrack, and from the
## 'groupAnnotation' parameter for an AnnotationTrack. The special value 'lowest' should work across all classes and provided the most
## detailed annotation ('id' for AnnotationTrack and 'exon' for GeneRegionTrack).

#' @describeIn AnnotationTrack-class return track item identifiers.
#' Depending on the setting of the optional argument lowest, these are either
#' the group identifiers or the individual item identifiers.
#' @export
setMethod("identifier", "AnnotationTrack", function(GdObject, type = .dpOrDefault(GdObject, "groupAnnotation", "group")) {
    if (is.null(type)) {
        type <- "group"
    }
    id <- switch(as.character(type),
        "group" = group(GdObject),
        "id" = .getAnn(GdObject, "id"),
        "feature" = feature(GdObject),
        "lowest" = .getAnn(GdObject, "id"),
        group(GdObject)
    )
    id[is.na(id)] <- "NA"
    return(id)
})

#' @describeIn AnnotationTrack-class Set the track item identifiers.
#' The replacement value has to be a character vector of appropriate length.
#' This always replaces the group-level identifiers, so essentially it is
#' similar to `groups<-`.
#' @export
setReplaceMethod("identifier", c("AnnotationTrack", "character"), function(GdObject, value) {
    type <- .dpOrDefault(GdObject, "groupAnnotation", "group")
    switch(as.character(type),
        "group" = group(GdObject) <- value,
        "id" = ranges(GdObject)$id <- value,
        "feature" = feature(GdObject) <- value,
        "lowest" = ranges(GdObject)$id <- value,
        group(GdObject) <- value
    )
    return(GdObject)
})
## Stacking ------------------------------------------------------------------

#' @describeIn AnnotationTrack-class Recompute the stacks based on the available
#' space and on the object's track items and stacking settings.
#' @export
setMethod("setStacks", "AnnotationTrack", function(GdObject, recomputeRanges = TRUE) {
    if (!.needsStacking(GdObject) || length(GdObject) == 0) {
        ## No stacks needed, so there is just a single bin
        bins <- rep(1, length(GdObject))
    } else {
        gp <- group(GdObject)
        if (!is.factor(gp)) {
            gp <- factor(gp, levels=unique(gp))
        }
        needsGrp <- any(duplicated(gp))
        lranges <- .dpOrDefault(GdObject, ".__groupRanges")
        gpt <- if (needsGrp) table(gp) else rep(1, length(GdObject))
        if (recomputeRanges || is.null(lranges) || length(lranges) != length(gpt)) {
            GdObject <- .computeGroupRange(GdObject)
            lranges <- .dpOrDefault(GdObject, ".__groupRanges")
        }
        uid <- if (.transcriptsAreCollapsed(GdObject)) {
            sprintf("uid%i", seq_along(identifier(GdObject)))
        } else {
            make.unique(identifier(GdObject, type = "lowest"))
        }

        bins <- rep(disjointBins(lranges), gpt)
        names(bins) <- if (needsGrp) unlist(split(uid, gp)) else uid
        bins <- bins[uid]
    }
    bins <- if (length(bins)) bins else 0
    GdObject@stacks <- bins
    return(GdObject)
})
## Consolidate ---------------------------------------------------------------

#' @describeIn AnnotationTrack-class Consolidate.
#' Determine whether there is group label annotation or not, and add this
#' information as the internal display parameter `.__hasAnno`. Precompute
#' the grouped ranges together with optional labels in order to determine
#' the correct plotting range later.
#' @keywords internal
setMethod("consolidateTrack", signature(GdObject = "AnnotationTrack"), function(GdObject, hasAxis = FALSE,
                                                                                hasTitle = .dpOrDefault(GdObject, "showTitle", TRUE),
                                                                                title.width = NULL, ...) {
    GdObject <- callNextMethod(GdObject, ...)
    if (length(GdObject)) {
        ## ids are shown if either set by the showId parameter, or through the use of the transcriptAnnotation or
        ## groupAnnotation parameters
        ids <- identifier(GdObject)
        sid <- .dpOrDefault(GdObject, "showId")
        ta <- if (is(GdObject, "GeneRegionTrack")) .dpOrDefault(GdObject, "transcriptAnnotation") else .dpOrDefault(GdObject, "groupAnnotation")
        sid <- if (is.null(sid)) !is.null(ta) && ta != "none" else sid
        hasAnno <- sid && !all(ids == "")
        ## element ids are shown if either the showFeatureId or showExonId parameters are TRUE, or if featureAnnotation
        ## or exonAnnotation is set
        sfid <- if (is(GdObject, "GeneRegionTrack")) .dpOrDefault(GdObject, "showExonId") else .dpOrDefault(GdObject, "showFeatureId")
        fa <- if (is(GdObject, "GeneRegionTrack")) .dpOrDefault(GdObject, "exonAnnotation") else .dpOrDefault(GdObject, "featureAnnotation")
        sfid <- if (is.null(sfid)) !is.null(fa) && fa != "none" else sfid
        displayPars(GdObject) <- list(".__hasAnno" = hasAnno, showId = sid, showFeatureId = sfid, showExonId = sfid)
        GdObject <- .computeGroupRange(GdObject, hasAxis = hasAxis, hasTitle = hasTitle, title.width = title.width)
    } else {
        displayPars(GdObject) <- list(".__hasAnno" = FALSE, showId = FALSE, showFeatureId = FALSE, showExonId = FALSE)
    }
    return(GdObject)
})

## Collapse  -----------------------------------------------------------------

## A slightly quicker function to compute overlaps between two GRanges objects
.myFindOverlaps <- function(gr1, gr2) {
    gr1 <- sort(gr1)
    gr2 <- sort(gr2)
    gr1 <- split(gr1, seqnames(gr1))
    gr2 <- split(gr2, seqnames(gr2))
    queryHits(findOverlaps(ranges(gr1), ranges(gr2)))
}
## Find all elements in a GRanges object 'grange' with distance smaller than 'minXDist' and merge them along with their additional
## metadata columns. 'elements' is a frequency table of items per group, and it is needed to figure out whether all items of a given
## group have been merged. 'GdObject' is the input track object from which certain information has to be extracted. The output of this
## function is a list with elements
##   o range: the updated GRanges object
##   o needsRestacking: logical flag indicating whether stacks have to be recomputed
##   o split: the original merged annotation (need to work on this, currently not used)
##   o merged: logical vector indicating which of the elements in 'range' constitute fully merged groups
.collapseAnnotation <- function(grange, minXDist, elements, GdObject, offset = 0) {
    needsRestacking <- TRUE
    annoSplit <- merged <- NULL
    anno <- as.data.frame(grange)
    for (i in colnames(anno)) {
        if (is.factor(anno[, i])) {
            anno[, i] <- as.character(anno[, i])
        }
    }
    cols <- c("strand", "density", "gdensity", "feature", "id", "start", "end", if (is(GdObject, "GeneRegionTrack")) {
        c("gene", "exon", "transcript", "symbol", "rank")
    } else {
        "group"
    })
    missing <- which(!cols %in% colnames(anno))
    for (i in missing) {
        anno[, cols[missing]] <- if (cols[i] == "density") 1 else NA
    }
    rRed <- if (length(grange) > 1) reduce(grange, min.gapwidth = minXDist, with.revmap = TRUE) else grange
    if (length(rRed) < length(grange)) {
        ## Some of the items have to be merged and we need to make sure that the additional annotation data that comes with it
        ## is processed in a sane way.
        needsRestacking <- TRUE
        mapping <- rep(seq_along(rRed$revmap), elementNROWS(rRed$revmap))
        ## We start by finding the items that have not been reduced
        identical <- mapping %in% which(table(mapping) == 1)
        newVals <- anno[identical, cols]
        ## Here we hijack the seqnames column to indicate whether the whole group has been merged
        if (nrow(newVals)) {
            newVals$seqnames <- elements[as.character(anno[identical, "seqnames"])] == 1
            newVals$gdensity <- ifelse(elements[as.character(anno[identical, "seqnames"])] == 1, 1, NA)
        }
        ## Now find out which original items have been merged
        grange <- grange[!identical]
        rRed <- rRed[-(mapping[identical])]
        index <- mapping[!identical]
        annoSplit <- split(anno[!identical, ], index)
        cid <- function(j) sprintf("[Cluster_%i]  ", j + offset)
        ## FIXME: We could speed this up by running it in C
        tmpVals <- lapply(seq_along(annoSplit), function(i) {
            x <- annoSplit[[i]]
            if (is(GdObject, "GeneRegionTrack")) {
                c(
                    strand = ifelse(length(unique(x[, "strand"])) == 1, as.character(x[1, "strand"]), "*"),
                    density = sum(as.integer(x[, "density"])),
                    gdensity = ifelse(is.na(head(x[, "gdensity"], 1)), 1, sum(as.integer(x[, "gdensity"]))),
                    feature = ifelse(length(unique(x[, "feature"])) == 1, as.character(x[1, "feature"]), "composite"),
                    id = ifelse(length(unique(x[, "id"])) == 1, as.character(x[1, "id"]), cid(i)),
                    start = min(x[, "start"]),
                    end = max(x[, "end"]),
                    gene = ifelse(length(unique(x[, "gene"])) == 1, as.character(x[1, "gene"]), cid(i)),
                    exon = ifelse(length(unique(x[, "exon"])) == 1, as.character(x[1, "exon"]), cid(i)),
                    transcript = ifelse(length(unique(x[, "transcript"])) == 1, as.character(x[1, "transcript"]), cid(i)),
                    symbol = ifelse(length(unique(x[, "symbol"])) == 1, as.character(x[1, "symbol"]), cid(i)),
                    rank = min(as.integer(x[, "rank"])), seqnames = as.vector(nrow(x) == elements[x[1, "seqnames"]])
                )
            } else {
                c(
                    strand = ifelse(length(unique(x[, "strand"])) == 1, as.character(x[1, "strand"]), "*"),
                    density = sum(as.integer(x[, "density"])),
                    gdensity = ifelse(is.na(head(x[, "gdensity"], 1)), 1, sum(as.integer(x[, "gdensity"]))),
                    feature = ifelse(length(unique(x[, "feature"])) == 1, as.character(x[1, "feature"]), "composite"),
                    id = ifelse(length(unique(x[, "id"])) == 1, as.character(x[1, "id"]), cid(i)),
                    start = min(x[, "start"]),
                    end = max(x[, "end"]),
                    group = ifelse(length(unique(x[, "group"])) == 1, as.character(x[1, "group"]), cid(i)),
                    seqnames = as.vector(nrow(x) == elements[x[1, "seqnames"]])
                )
            }
        })
        newVals <- rbind(newVals, as.data.frame(do.call(rbind, tmpVals), stringsAsFactors = FALSE))
        merged <- as.logical(newVals$seqnames)
        grange <- GRanges(
            seqnames = chromosome(GdObject), strand = newVals[, "strand"],
            ranges = IRanges(start = as.integer(newVals[, "start"]), end = as.integer(newVals[, "end"]))
        )
        cnMatch <- match(c(colnames(values(GdObject)), "gdensity"), colnames(newVals))
        mcols(grange) <-
            if (any(is.na(cnMatch))) newVals[, setdiff(colnames(newVals), c("strand", "start", "end", "seqnames"))] else newVals[, cnMatch]
    } else {
        grange2 <- GRanges(seqnames = chromosome(GdObject), strand = strand(grange), ranges = ranges(grange))
        mcols(grange2) <- mcols(grange)
        grange <- grange2
    }
    return(list(range = grange, needsRestacking = needsRestacking, split = annoSplit, merged = merged, offset = length(annoSplit)))
}

## For AnnotationTracks we need to collapse the all regions along with the additional annotation.
## For GeneRegionTracks we essentially need to do the same thing as for AnnotationTracks, however the additional annotation columns
## are quite different. We do this in multiple turn with increasing levels of complexity:
##    1.) merge all individual items within a group that can no longer be separated
##    2.) merge overlapping groups with just a single remaining item (optional, if mergeGroups==TRUE)

#' @describeIn AnnotationTrack-class preprocess the track before plotting.
#' This will collapse overlapping track items based on the available resolution
#' and increase the width and height of all track objects to a minimum value
#' to avoid rendering issues. See collapsing for details.
#' @keywords internal
setMethod(
    "collapseTrack", signature(GdObject = "AnnotationTrack"),
    function(GdObject, diff = .pxResolution(coord = "x"), xrange) {
        ## We first add the original unmodified GdObject as a display parameter to be able to reference back if we ever need to
        displayPars(GdObject) <- list(".__OriginalGdObject" = .deepCopyPars(GdObject))
        collapse <- .dpOrDefault(GdObject, "collapse", TRUE)
        min.width <- .dpOrDefault(GdObject, "min.width", 2)
        min.distance <- .dpOrDefault(GdObject, "min.distance", 2)
        minXDist <- max(0, ceiling(min.distance * diff))
        r <- ranges(GdObject)
        ## Compute native coordinate equivalent to 1 pixel and resize
        rNew <- .resize(r, min.width, diff)
        needsRestacking <- any(r != rNew)
        r <- rNew
        ## Collapse all items within a group to a single meta-item (if collapseTranscripts is set)
        if (is(GdObject, "GeneRegionTrack")) {
            ctrans <- .dpOrDefault(GdObject, "collapseTranscripts", FALSE)
            if (is.logical(ctrans) && ctrans) {
                ctrans <- "gene"
            }
            switch(ctrans,
                "gene" = {
                    newVals <- unlist(endoapply(split(values(r), paste(gene(GdObject), strand(GdObject))), head, 1))
                    newVals$exon <- NA
                    newVals$feature <- "merged"
                    newVals$transcript <- newVals$gene
                    r <- unlist(range(split(r, gene(GdObject))))
                    mcols(r) <- newVals
                    GdObject@range <- r
                },
                "longest" = {
                    r <- unlist(endoapply(split(r, gene(GdObject)), function(x) {
                        xs <- split(x, x$transcript)
                        xs[[which.max(vapply(xs, function(y) abs(diff(as.numeric(as.data.frame(ranges(range(y)))[, c("start", "end")]))), FUN.VALUE = numeric(1L)))]]
                    }))
                    GdObject@range <- r
                },
                "shortest" = {
                    r <- unlist(endoapply(split(r, gene(GdObject)), function(x) {
                        xs <- split(x, x$transcript)
                        xs[[which.min(vapply(xs, function(y) abs(diff(as.numeric(as.data.frame(ranges(range(y)))[, c("start", "end")]))), FUN.VALUE = numeric(1L)))]]
                    }))
                    GdObject@range <- r
                },
                "meta" = {
                    newVals <- unlist(endoapply(split(values(r), paste(gene(GdObject), strand(GdObject))), head, 1))
                    newVals$feature <- "merged"
                    newVals$transcript <- newVals$gene
                    rtmp <- reduce(split(r, paste(gene(GdObject), strand(GdObject))))
                    newVals <- newVals[rep(seq_along(rtmp), elementNROWS(rtmp)), ]
                    newVals$exon <- paste("merged_exon_", unlist(lapply(elementNROWS(rtmp), function(x) seq(1, x)), use.names = FALSE), sep = "")
                    r <- unlist(rtmp)
                    mcols(r) <- newVals
                    GdObject@range <- r
                }
            )
        }
        ## Collapse overlapping ranges (less than minXDist space between them) and process the annotation data
        if (collapse) {
            ## Merge all items in those groups for which no individual items can be separated
            elements <- table(group(GdObject))
            rr <- GRanges(seqnames = as.character(group(GdObject)), ranges = IRanges(start = start(r), end = end(r)), strand = strand(r))
            mcols(rr) <- mcols(r)
            rr <- sort(unique(rr))
            mergedAnn <- .collapseAnnotation(rr, minXDist, elements, GdObject)
            needsRestacking <- needsRestacking || mergedAnn$needsRestacking
            ## Now we take a look whether there are any groups that could be merged (if mergeGroups is TRUE)
            if (.dpOrDefault(GdObject, "mergeGroups", FALSE) && any(mergedAnn$merged)) {
                rr <- sort(mergedAnn$range[mergedAnn$merged])
                strand(rr) <- "*"
                mergedAnn2 <- .collapseAnnotation(rr, minXDist, elements, GdObject, mergedAnn$offset)
                needsRestacking <- needsRestacking || mergedAnn2$needsRestacking
                mergedAnn$range <- c(mergedAnn$range[!mergedAnn$merged], mergedAnn2$range)
            }
            r <- mergedAnn$range
        }
        ## Reconstuct the track object and return
        GdObject@range <- r
        ## if(needsRestacking)
        GdObject <- setStacks(GdObject)
        return(GdObject)
    }
)

## Subset --------------------------------------------------------------------

## In order to keep the grouping information for track regions in the clipped areas we have to
## keep all group elements that overlap with the range. We still want to record the requested
## ranges in the internal '.__plottingRange' display parameter.

#' @describeIn AnnotationTrack-class subset a `AnnotationTrack` by coordinates
#' and sort if necessary.
#' @export
setMethod("subset", signature(x = "AnnotationTrack"), function(x, from = NULL, to = NULL, sort = FALSE, stacks = FALSE, use.defaults = TRUE, ...) {
    ## Subset to a single chromosome first
    lx <- length(x)
    csel <- seqnames(x) != chromosome(x)
    if (any(csel)) {
        x <- x[!csel]
    }
    if (length(x)) {
        ## Nothing to do if everything is within the range
        granges <- unlist(range(split(ranges(x), group(x))))
        ranges <- if (use.defaults) {
            .defaultRange(x, from = from, to = to)
        } else {
            c(
                from = ifelse(is.null(from), min(start(granges)) - 1, from),
                to = ifelse(is.null(to), max(end(granges)) + 1, to)
            )
        }
        if (!(any(end(x) < ranges["from"] | start(x) > ranges["to"]))) {
            if (stacks) {
                x <- setStacks(x)
            }
            return(x)
        }
        ## Now remove everything except for the overlapping groups by first subselecting all groups in the range...
        gsel <- names(granges)[subjectHits(findOverlaps(GRanges(seqnames = chromosome(x), ranges = IRanges(min(ranges), max(ranges))), granges))]
        x <- x[group(x) %in% gsel]
        if (sort) {
            x <- x[order(range(x)), ]
        }
        if (stacks) {
            x <- setStacks(x)
        }
        displayPars(x) <- list(".__plottingRange" = ranges)
    }
    if (length(x) != lx) {
        x <- .computeGroupRange(x)
    }
    return(x)
})

## ReferenceDataTracks need to stream the data from file and then pass the results on to the next method

#' @describeIn AnnotationTrack-class subset a `ReferenceAnnotationTrack` by
#' coordinates and sort if necessary.
#' @export
setMethod("subset", signature(x = "ReferenceAnnotationTrack"), function(x, from, to, chromosome, ...) {
    ## We only need to reach out into the referenced file once if the range is already contained in the object
    if (missing(from) || is.null(from) || missing(to) || is.null(to)) {
        stop("Need both start and end location to subset a ReferenceAnnotationTrack")
    }
    if (missing(chromosome) || is.null(chromosome)) {
        chromosome <- Gviz::chromosome(x)
    }
    subRegion <- GRanges(seqnames = chromosome[1], ranges = IRanges(start = from, end = to))
    if (length(ranges(x)) == 0 || all(overlapsAny(ranges(x), subRegion))) {
        cMap <- .resolveColMapping(x@stream(x@reference, subRegion), x@args, x@mapping)
        x@range <- .buildRange(cMap$data, args = cMap$args, defaults = x@defaults, trackType = "AnnotationTrack")
        chromosome(x) <- chromosome[1]
    }
    return(callNextMethod(x = x, from = from, to = to, drop = FALSE, ...))
})

## Position ------------------------------------------------------------------
## DrawGrid ------------------------------------------------------------------
setMethod("drawGrid", signature(GdObject = "AnnotationTrack"), function(GdObject, from, to) {
    if (.dpOrDefault(GdObject, "grid", FALSE)) {
        pushViewport(dataViewport(xData = c(from, to), extension = c(0, 0), yData = 0:1, clip = TRUE))
        panel.grid(
            h = 0, v = .dpOrDefault(GdObject, "v", -1),
            col = .dpOrDefault(GdObject, "col.grid", "#e6e6e6"), lty = .dpOrDefault(GdObject, "lty.grid", 1),
            lwd = .dpOrDefault(GdObject, "lwd.grid", 1)
        )
        popViewport(1)
    }
})
## DrawAxis ------------------------------------------------------------------
## DrawGD --------------------------------------------------------------------

## Draw gene models as found in AnnotationTracks or GeneRegionTracks
##
## Calculate all coordinates and values for the individual stacks first, and append to the
## input list 'toDraw'. This allows us to plot the whole track at once, making use of grid's
## vectorization. Without this tweak, every stack would have to be drawn individually, which
## can be painfully slow...
## Because all of the y coordinates are calculated in the coordinate system of the current
## viewport for the respective stack we have to add offset values.

## Compute the coordinates and colors for all track items (e.g. exons in a GeneRegionTrack)
.boxes <- function(GdObject, offsets) {
    ylim <- c(0, 1)
    h <- diff(ylim)
    middle <- mean(ylim)
    sh <- max(0, min(h, .dpOrDefault(GdObject, "stackHeight", 0.75)))
    space <- (h - (h * sh)) / 2
    shift <- switch(.dpOrDefault(GdObject, "stackJust", "middle"),
        "top" = space,
        "bottom" = -space,
        0
    )
    if (inherits(GdObject, "GeneRegionTrack")) {
        thinBox <- .dpOrDefault(GdObject, "thinBoxFeature", .THIN_BOX_FEATURES)
        space <- ifelse(feature(GdObject) %in% thinBox, space + ((middle -
            space) / 2), space)
    }
    shape <- .dpOrDefault(GdObject, "shape", "arrow")
    color <- .getBiotypeColor(GdObject)
    id <- identifier(GdObject, type = .dpOrDefault(
        GdObject, ifelse(is(GdObject, "GeneRegionTrack"), "exonAnnotation", "featureAnnotation"),
        "lowest"
    ))
    sel <- grepl("\\[Cluster_[0-9]*\\]", id)
    id[sel] <- sprintf(
        "%i merged\n%s", as.integer(.getAnn(GdObject, "density")[sel]),
        ifelse(class(GdObject) %in% c("AnnotationTrack", "DetailsAnnotationTrack"), "features", "exons")
    )
    yy1 <- ylim[1] + space + offsets + shift
    yy2 <- ylim[2] - space + offsets + shift
    boxes <- data.frame(
        cx1 = start(GdObject), cy1 = yy1, cx2 = start(GdObject) + width(GdObject), cy2 = yy2,
        fill = color, strand = strand(GdObject), text = id, textX = start(GdObject) + (width(GdObject) / 2), textY = middle + offsets,
        .getImageMap(cbind(start(GdObject), yy1, end(GdObject), yy2)),
        start = start(GdObject), end = end(GdObject), values(GdObject), exonId = id, origExonId = identifier(GdObject, type = "lowest"),
        stringsAsFactors = FALSE
    )
    rownames(boxes) <- if (.transcriptsAreCollapsed(GdObject)) {
        sprintf("uid%i", seq_along(identifier(GdObject)))
    } else {
        make.unique(identifier(GdObject, type = "lowest"))
    }
    return(boxes)
}

## Compute the coordinates for the bars connecting grouped items and the group labels

#' @importFrom grDevices col2rgb hsv rgb2hsv
.barsAndLabels <- function(GdObject) {
    bins <- stacks(GdObject)
    stacks <- max(bins)
    res <- .pxResolution(coord = "x")
    gp <- group(GdObject)
    if (!is.factor(gp)) {
        gp <- factor(gp, levels=unique(gp))
    }
    grpSplit <- split(range(GdObject), gp)
    grpRanges <- unlist(range(grpSplit))
    needBar <- vapply(grpSplit, length, FUN.VALUE = numeric(1L)) > 1 & width(grpRanges) > res
    ## If we draw the bar from start to end of the range we sometimes see little overlaps that extend beyond the first or last item.
    ## In order to fix this, we just subtract the equivalent of min.width pixels from both ends of each group range
    min.swidth <- res * .dpOrDefault(GdObject, "min.width", 2)
    nstart <- start(grpRanges[needBar]) + min.swidth
    nend <- end(grpRanges[needBar]) - min.swidth
    sel <- (nend - nstart) > 0
    start(grpRanges[needBar][sel]) <- nstart[sel]
    end(grpRanges[needBar][sel]) <- nend[sel]
    strand <- vapply(split(strand(GdObject), gp), function(x) {
        tmp <- unique(x)
        if (length(tmp) > 1) "*" else tmp
    }, FUN.VALUE = character(1L))
    yloc <- vapply(split((stacks - bins) + 1, gp), function(x) unique(x), FUN.VALUE = numeric(1L)) + 0.5
    color <- if (length(grep("__desatCol", values(GdObject)$feature[1]))) {
        .dpOrDefault(GdObject, "fill", .DEFAULT_FILL_COL)
    } else {
        vapply(split(.getBiotypeColor(GdObject), gp), head, 1, FUN.VALUE = character(1L))
    }
    bars <- data.frame(
        sx1 = start(grpRanges)[needBar], sx2 = end(grpRanges)[needBar], y = yloc[needBar], strand = strand[needBar],
        col = color[needBar], stringsAsFactors = FALSE
    )
    labs <- .dpOrDefault(GdObject, ".__groupLabels")[names(grpRanges)]
    if (!is.null(labs)) {
        lsel <- grepl("\\[Cluster_[0-9]*\\]", labs)
        if (any(lsel)) {
            gdens <- as.integer(vapply(split(.getAnn(GdObject, "gdensity"), gp), head, 1, FUN.VALUE = character(1L)))
            labs[lsel] <- sprintf(
                "%i merged %s  ", gdens[lsel],
                ifelse(class(GdObject) %in% c("AnnotationTrack", "DetailsAnnotationTrack"), "groups", "transcript models")
            )
        }
        just <- .dpOrDefault(GdObject, "just.group", "left")
        rev <- .dpOrDefault(GdObject, "reverseStrand", FALSE)
        sizes <- .dpOrDefault(GdObject, ".__groupLabelWidths")[names(grpRanges), , drop = FALSE]
        pr <- .dpOrDefault(GdObject, ".__plottingRange", data.frame(from = min(start(GdObject)), to = max(end(GdObject))))
        ## grpRangesCut <- restrict(grpRanges, start=as.integer(pr["from"]), end=as.integer(pr["to"]))
        switch(just,
            "left" = {
                cx <- if (!rev) start(grpRanges) - sizes$after else end(grpRanges) + sizes$after
                cy <- yloc
                algn <- c("right", "center")
            },
            "right" = {
                cx <- if (!rev) end(grpRanges) + sizes$after else start(grpRanges) - sizes$after
                cy <- yloc
                algn <- c("left", "center")
            },
            "above" = {
                ## cx <- start(grpRangesCut) + width(grpRangesCut)/2
                ## indx <- which(seq_len(length(grpRanges)) %in% queryHits(findOverlaps(grpRanges, grpRangesCut)))
                ## cx <- cx[indx] # needs revision
                ## cy <- yloc[indx] + 0.5
                ## labs <- labs[indx]
                cx <- start(grpRanges) + width(grpRanges) / 2
                cy <- yloc + 0.5
                algn <- c("center", "top")
            },
            "below" = {
                ## cx <- start(grpRangesCut) + width(grpRangesCut)/2
                ## indx <- which(seq_len(length(grpRanges)) %in% queryHits(findOverlaps(grpRanges, grpRangesCut)))
                ## cx <- cx[indx] # needs revision
                ## cy <- yloc[indx] - 0.5
                ## labs <- labs[indx]
                cx <- start(grpRanges) + width(grpRanges) / 2
                cy <- yloc - 0.5
                algn <- c("center", "bottom")
            },
            stop(sprintf("Unknown label justification '%s'", just))
        )

        labels <- data.frame(txt = labs, x = cx, y = cy, stringsAsFactors = FALSE)
    } else {
        labels <- algn <- NA
    }
    return(list(bars = bars, labels = labels, align = algn))
}

#' @describeIn AnnotationTrack-class plot the object to a graphics device.
#' The return value of this method is the input object, potentially updated
#' during the plotting operation. Internally, there are two modes in which the
#' method can be called. Either in 'prepare' mode, in which case no plotting is
#' done but the object is preprocessed based on the available space, or in
#' 'plotting' mode, in which case the actual graphical output is created.
#' Since subsetting of the object can be potentially costly, this can be
#' switched off in case subsetting has already been performed before or
#' is not necessary.
#'
#' @export
setMethod("drawGD", signature("AnnotationTrack"), function(GdObject, minBase, maxBase, prepare = FALSE, subset = TRUE, ...) {
    debug <- .dpOrDefault(GdObject, "debug", FALSE)
    if ((is.logical(debug) && debug) || debug == "prepare") {
        browser()
    }
    imageMap(GdObject) <- NULL
    if (!length(GdObject)) {
        return(invisible(GdObject))
    }
    ## In prepare mode we need to make sure that the stacking information is updated from the optional display parameter (by calling
    ## the StackedTrack drawGD method) and also perform the collapsing of track items which could potentially lead to re-stacking.
    if (prepare) {
        GdObject <- callNextMethod(GdObject, ...)
        bins <- stacks(GdObject)
        stacks <- max(bins)
        ## We need to collapse the track object based on the current screen resolution (note that this may trigger re-stacking)
        pushViewport(dataViewport(xData = c(minBase, maxBase), extension = 0, yscale = c(1, stacks + 1), clip = TRUE))
        GdObject <- collapseTrack(GdObject, diff = .pxResolution(coord = "x"), xrange = c(minBase, maxBase))
        popViewport(1)
        return(invisible(GdObject))
    }
    if ((is.logical(debug) && debug) || debug == "draw") {
        browser()
    }
    ## If there are too many stacks for the available device resolution we cast an error, otherwise we set the viewport
    bins <- stacks(GdObject)
    stacks <- max(bins)
    rev <- .dpOrDefault(GdObject, "reverseStrand", FALSE)
    xscale <- if (!rev) c(minBase, maxBase) else c(maxBase, minBase)
    yscale <- if (!.dpOrDefault(GdObject, "reverseStacking", FALSE)) c(1, stacks + 1) else c(stacks + 1, 1)
    pushViewport(dataViewport(xscale = xscale, extension = 0, yscale = yscale, clip = TRUE))
    res <- .pxResolution(coord = "x")
    curVp <- vpLocation()
    if (curVp$size["height"] / stacks < .dpOrDefault(GdObject, "min.height", 3)) {
        stop("Too many stacks to draw. Either increase the device size or limit the drawing to a smaller region.")
    }
    ## We adjust the color saturation to indicate overplotting if necessary
    if (.dpOrDefault(GdObject, "showOverplotting", FALSE)) {
        dens <- as.numeric(values(GdObject)$density)
        if (length(unique(dens)) != 1) {
            minSat <- max(0.25, 1 / max(dens))
            minDens <- min(dens)
            rDens <- diff(range(dens))
            saturation <- minSat + ((dens - minDens) / rDens / (1 / (1 - minSat)))
            bc <- unique(.getBiotypeColor(GdObject))
            baseCol <- rgb2hsv(col2rgb(bc))
            desatCols <- unlist(lapply(saturation, function(x) hsv(baseCol[1, ], x, baseCol[3, ])))
            names(desatCols) <- paste(unique(feature(GdObject)), rep(dens, each = length(bc)), sep = "_")
            feature(GdObject) <- paste(feature(GdObject), dens, sep = "_")
            desatCols <- desatCols[unique(names(desatCols))]
            displayPars(GdObject) <- as.list(desatCols)
        }
    }
    ## Now we can pre-compute all the coordinates and settings for the elements to be drawn, ...
    box <- .boxes(GdObject, (stacks - bins) + 1)
    barsAndLab <- .barsAndLabels(GdObject)
    bar <- barsAndLab$bars
    bartext <- barsAndLab$labels
    ## ... get all the necessary display parameters
    shape <- .dpOrDefault(GdObject, "shape", "arrow")
    col.line <- .dpOrDefault(GdObject, "col.line")[1]
    border <- .dpOrDefault(GdObject, "col")[1]
    if (is.null(border)) {
        border <- ifelse(is(GdObject, "GeneRegionTrack"), NA, "transparent")
    }
    lwd <- .dpOrDefault(GdObject, "lwd", 2)
    lty <- .dpOrDefault(GdObject, "lty", 1)
    alpha <- .dpOrDefault(GdObject, "alpha", 1)
    rotation <- .dpOrDefaultFont(GdObject, "rotation", "item", 0)
    rotation.group <- .dpOrDefaultFont(GdObject, "rotation", "group", 0)
    just <- .dpOrDefault(GdObject, "just.group", "left")
    ## ... and finally draw whatever is needed
    if (nrow(box) > 0) {
        ## If we want to place a label on top or below the ranges we need to know how much size that will take up
        ## and adjust the plotting shapes accordingly
        drawLabel <- .dpOrDefault(GdObject, ".__hasAnno", FALSE) && !is.null(bartext) && !anyNA(bartext) &&
            nrow(bartext) > 0 && stacking(GdObject) != "dense"
        bs <- as.vector((stacks - bins) + 1)
        if (drawLabel && just %in% c("above", "below")) {
            labelHeights <- max(.getStringDims(GdObject, bartext$txt, subtype = "group")$height)
            labelSpace <- as.numeric(convertHeight(unit(2, "points"), "native"))
            avSpace <- (1 - max(box$cy2 - box$cy1)) / 2
            vadjust <- max(0, (labelHeights + (2 * labelSpace)) - avSpace)
            bh <- 1 - (avSpace * 2)
            bhNew <- bh - vadjust
            sfac <- bhNew / bh
            if (just == "above") {
                bar$y <- bar$y - (vadjust / 2)
                box$textY <- box$textY - (vadjust / 2)
                bartext$y <- bartext$y - labelSpace
                box$cy1 <- bs + ((box$cy1 %% 1 - avSpace) * sfac) + avSpace
                box$cy2 <- bs + ((box$cy2 %% 1 - avSpace) * sfac) + avSpace
            } else {
                bar$y <- bar$y + (vadjust / 2)
                box$textY <- box$textY + (vadjust / 2)
                bartext$y <- bartext$y + labelSpace
                box$cy1 <- bs + ((box$cy1 %% 1 - avSpace) * sfac) + avSpace + vadjust
                box$cy2 <- bs + ((box$cy2 %% 1 - avSpace) * sfac) + avSpace + vadjust
            }
        }
        ## Plotting of the (arrow)bar
        if (nrow(bar) > 0) {
            .arrowBar(bar$sx1, bar$sx2,
                y = bar$y, bar$strand, box[, seq_len(4), drop = FALSE],
                W = .dpOrDefault(GdObject, "arrowFeatherWidth", 3), D = .dpOrDefault(GdObject, "arrowFeatherDistance", 20),
                col = if (is.null(col.line)) bar$col else rep(col.line, length(bar$col)), lwd = lwd, lty = lty,
                alpha = alpha, barOnly = (!"smallArrow" %in% .dpOrDefault(GdObject, "shape", "box") || stacking(GdObject) == "dense"),
                diff = res, min.height = .dpOrDefault(GdObject, "min.height", 3)
            )
        }
        ## Plotting of the boxes
        box$col <- if (is.na(border)) box$fill else border
        if ("box" %in% shape || ("smallArrow" %in% shape && !("arrow" %in% shape || "ellipse" %in% shape))) {
            .filledBoxes(box, lwd = lwd, lty = lty, alpha = alpha)
        }
        ## Plotting of the ellipses
        if ("ellipse" %in% shape) {
            ellCoords <- .box2Ellipse(box)
            grid.polygon(
                x = ellCoords$x1, y = ellCoords$y1, id = ellCoords$id,
                gp = gpar(col = if (is.na(border)) box$fill else border, fill = box$fill, lwd = lwd, lty = lty, alpha = alpha),
                default.units = "native"
            )
        }
        ## Plotting of the filled arrows
        if ("arrow" %in% shape && !"box" %in% shape) {
            .filledArrow(box, lwd = lwd, lty = lty, alpha = alpha, min.width = 4 * res, max.width = .dpOrDefault(GdObject, "arrowHeadMaxWidth", 40) * res)
        }
        ## Plotting of the filled arrows with fixed head size
        if ("fixedArrow" %in% shape && !"box" %in% shape) {
            .filledArrow(box,
                lwd = lwd, lty = lty, alpha = alpha, min.width = 4 * res, absoluteWidth = TRUE,
                W = .dpOrDefault(GdObject, "arrowHeadWidth", 30) * res
            )
        }
        ## Plotting of the item labels
        if (.dpOrDefault(GdObject, "showFeatureId", FALSE)) {
            # grid.text(str2expression(box$text), box$textX, box$textY,
            grid.text(box$text, box$textX, box$textY,
                rot = rotation, gp = .fontGp(GdObject, subtype = "item"),
                default.units = "native", just = c("center", "center")
            )
        }
        ## Plotting of the group labels
        if (drawLabel) {
            grid.text(bartext$txt, bartext$x, bartext$y,
                rot = rotation.group, gp = .fontGp(GdObject, subtype = "group"),
                default.units = "native", just = barsAndLab$align
            )
        }
    }
    popViewport(1)
    ## Finally we set up the image map
    ## FIXME: we may want to record the merging information here
    im <- if (!is.null(box)) {
        coords <- as.matrix(box[, c("x1", "y1", "x2", "y2"), drop = FALSE])
        restCols <- setdiff(colnames(box), c("x1", "x2", "y1", "y2", "cx1", "cx2", "cy1", "cy2", "textX", "textY"))
        tags <- lapply(restCols, function(x) {
            tmp <- as.character(box[, x])
            names(tmp) <- rownames(coords)
            tmp
        })
        names(tags) <- restCols
        tags$title <- identifier(GdObject)
        ImageMap(coords = coords, tags = tags)
    } else {
        NULL
    }
    imageMap(GdObject) <- im
    return(invisible(GdObject))
})

## DrawGD DetailsAnnotationTrack ---------------------------------------------

## Create a data.frame with the distinct details function arguments (like start, end, ...)
.buildArgsDf <- function(GdObject) {
    groupDetails <- .dpOrDefault(GdObject, "groupDetails", FALSE)
    rr <- if (groupDetails) unlist(range(split(ranges(GdObject), group(GdObject)))) else ranges(GdObject)
    args <- data.frame(
        start = as.integer(start(rr)), end = as.integer(end(rr)), strand = as.character(strand(rr)),
        chromosome = as.character(seqnames(rr)),
        identifier = as.character(if (groupDetails) names(rr) else identifier(GdObject, type = "lowest")),
        stringsAsFactors = FALSE
    )
    return(args)
}

#' @describeIn AnnotationTrack-class plot the object to a graphics device.
#' The return value of this method is the input object, potentially updated
#' during the plotting operation. Internally, there are two modes in which the
#' method can be called. Either in 'prepare' mode, in which case no plotting is
#' done but the object is preprocessed based on the available space, or in
#' 'plotting' mode, in which case the actual graphical output is created.
#' Since subsetting of the object can be potentially costly, this can be
#' switched off in case subsetting has already been performed before or
#' is not necessary.
#'
#' @export
setMethod(
    "drawGD", signature("DetailsAnnotationTrack"),
    function(GdObject, minBase, maxBase, prepare = FALSE, ...) {
        debug <- .dpOrDefault(GdObject, "debug", FALSE)
        if ((is.logical(debug) && debug) || debug == "prepare") {
            browser()
        }
        adf <- .buildArgsDf(GdObject)
        args <- .dpOrDefault(GdObject, "detailsFunArgs", fromPrototype = TRUE)
        groupDetails <- .dpOrDefault(GdObject, "groupDetails", FALSE)
        if (prepare) {
            GdObject <- callNextMethod()
            GdObject <- GdObject[order(start(GdObject))]
            indices <- if (groupDetails) seq_len(length(unique(group(GdObject)))) else seq_len(length(GdObject))
            xscale <- if (!.dpOrDefault(GdObject, "reverseStrand", FALSE)) c(minBase, maxBase) else c(maxBase, minBase)
            pushViewport(viewport(xscale = xscale))
            select <- lapply(indices, function(i) {
                warn <- FALSE
                iargs <- as.list(adf[i, ])
                iargs$index <- i
                iargs$GdObject <- GdObject
                iargs$GdObject.original <- .dpOrDefault(GdObject, ".__OriginalGdObject", GdObject)
                args <- c(args[setdiff(names(args), names(iargs))], iargs)
                res <- do.call(GdObject@selectFun, args)
                if (length(res) != 1 || !is.logical(res) || is.na(res)) {
                    res <- TRUE
                    warn <- TRUE
                }
                c(res = res, warn = warn)
            })
            select <- do.call(rbind, select)
            if (any(select[, "warn"])) {
                warning("The result of function 'selectFun' has to be a single logical value. Forcing the value to 'TRUE'")
            }
            select <- select[, "res"]
            popViewport(1)
            displayPars(GdObject) <- list(".__select" = select)
            return(invisible(GdObject))
        }
        if ((is.logical(debug) && debug) || debug == "draw") {
            browser()
        }
        n <- length(GdObject)
        col <- rep(.dpOrDefault(GdObject, "detailsConnector.col", fromPrototype = TRUE), n)[seq_len(n)]
        lty <- rep(.dpOrDefault(GdObject, "detailsConnector.lty", fromPrototype = TRUE), n)[seq_len(n)]
        lwd <- rep(.dpOrDefault(GdObject, "detailsConnector.lwd", fromPrototype = TRUE), n)[seq_len(n)]
        pch <- rep(.dpOrDefault(GdObject, "detailsConnector.pch", fromPrototype = TRUE), n)[seq_len(n)]
        cex <- rep(.dpOrDefault(GdObject, "detailsConnector.cex", fromPrototype = TRUE), n)[seq_len(n)]
        border.lty <- rep(.dpOrDefault(GdObject, "detailsBorder.lty", fromPrototype = TRUE), n)[seq_len(n)]
        border.lwd <- rep(.dpOrDefault(GdObject, "detailsBorder.lwd", fromPrototype = TRUE), n)[seq_len(n)]
        border.col <- rep(.dpOrDefault(GdObject, "detailsBorder.col", fromPrototype = TRUE), n)[seq_len(n)]
        border.fill <- rep(.dpOrDefault(GdObject, "detailsBorder.fill", fromPrototype = TRUE), n)[seq_len(n)]
        minwidth <- .dpOrDefault(GdObject, "details.minWidth", fromPrototype = TRUE)
        size <- .dpOrDefault(GdObject, "details.size", fromPrototype = TRUE)
        xyratio <- .dpOrDefault(GdObject, "details.ratio", fromPrototype = TRUE)
        if (0 >= size || size > 1) {
            warning("details.size must be >0 and <1 - reset to 0.5")
            size <- 0.5
        }
        selection <- .dpOrDefault(GdObject, ".__select", rep(TRUE, length(GdObject)))
        len <- sum(selection)
        bins <- if (!groupDetails) stacks(GdObject) else unlist(lapply(split(stacks(GdObject), group(GdObject)), unique))
        stacks <- max(bins)
        if (len > 0) {
            if (((maxBase - minBase) / len) / .pxResolution(coord = "x") < minwidth) {
                warning("too much detail for available space (plot fewer annotation or increase details.minWidth)!")
                popViewport(1)
                GdObject <- callNextMethod()
                return(GdObject)
            }
            rr <- if (groupDetails) unlist(range(split(ranges(GdObject), group(GdObject)))) else ranges(GdObject)
            xloc1 <- (end(rr) - start(rr)) / 2 + start(rr)
            yloc1 <- (stacks - (bins - 0.5) + 1)
            xloc2 <- ((1 / len * seq_len(len)) - 1 / len + (1 / len * 0.5))
            yloc2 <- rep(1, len)
            ## draw details plots (via user supplied function 'fun')
            pushViewport(viewport(height = size, y = 1 - size, just = c(0.5, 0)))
            w <- 1
            v <- 0
            vpl <- vpLocation()
            r <- vpl$size["width"] / len / vpl$size["height"]
            if (r > xyratio) {
                w <- xyratio / r
                v <- ((1 / len) - (1 / len * w)) / 2
            }
            indices <- if (groupDetails) seq_len(length(unique(group(GdObject)))) else seq_len(length(GdObject))
            j <- 1
            pres <- list()
            hasError <- FALSE
            for (i in indices[selection]) {
                pushViewport(viewport(width = 1 / len * w, x = ((1 / len * j) - 1 / len) + (v), just = c(0, 0.5)))
                grid.rect(gp = gpar(col = border.col[i], lwd = border.lwd[i], lty = border.lty[i], fill = border.fill[i]))
                iargs <- as.list(adf[i, ])
                iargs$index <- i
                iargs$GdObject <- GdObject
                iargs$GdObject.original <- .dpOrDefault(GdObject, ".__OriginalGdObject", GdObject)
                args <- c(args[setdiff(names(args), names(iargs))], iargs)
                pres[[as.character(j)]] <- try(do.call(GdObject@fun, args), silent = TRUE)
                if (!is.null(pres) && is(pres[[as.character(j)]], "try-error")) {
                    hasError <- TRUE
                    grid.segments(x0 = c(0.1, 0.1), x1 = c(0.9, 0.9), y0 = c(0.9, 0.1), y1 = c(0.1, 0.9), gp = gpar(col = "red", lwd = 3))
                }
                popViewport(1)
                j <- j + 1
            }
            if (hasError) {
                warning("There have been errors in the detail plotting function:\n", paste(pres, collapse = "\n"))
            }
            popViewport(1)
            ## plot AnnotationTrack and connectors to details
            xscale <- if (!.dpOrDefault(GdObject, "reverseStrand", FALSE)) c(minBase, maxBase) else c(maxBase, minBase)
            pushViewport(viewport(
                xscale = xscale,
                yscale = c(1, stacks + 1), clip = FALSE,
                height = 1 - size, y = 0, just = c(.5, 0)
            ))
            GdObject <- callNextMethod()
            grid.segments(
                x0 = unit(xloc1[selection], "native"), x1 = xloc2, y0 = unit(yloc1[selection], "native"),
                y1 = yloc2, gp = gpar(col = col, lwd = lwd, lty = lty, cex = cex)
            )
            grid.points(x = unit(xloc2, "npc"), y = unit(yloc2, "npc"), gp = gpar(col = col, cex = cex), pch = pch)
            grid.points(x = unit(xloc1[selection], "native"), y = unit(yloc1[selection], "native"), gp = gpar(col = col, cex = cex), pch = pch)
            popViewport(1)
        } else {
            GdObject <- callNextMethod()
        }
        return(invisible(GdObject))
    }
)

## SetAs ---------------------------------------------------------------------

#' @importFrom rtracklayer GenomicData
setAs(
    "AnnotationTrack", "UCSCData",
    function(from, to) {
        ranges <- range(from)
        dcolor <- as.integer(col2rgb(.dpOrDefault(from, "col")))
        line <- new("BasicTrackLine",
            name = names(from),
            description = names(from),
            visibility = stacking(from), color = dcolor, itemRgb = TRUE
        )
        vals <- values(from)
        color <- .getBiotypeColor(from)
        strand <- as.character(strand(from))
        strand[strand == "*"] <- "+"
        new("UCSCData", rtracklayer::GenomicData(ranges,
            chrom = chromosome(from),
            id = gsub(" ", "_", vals$id),
            name = gsub(" ", "_", as.character(vals$id)), itemRgb = color,
            strand = strand
        ),
        trackLine = line
        )
    }
)

setAs("GRanges", "AnnotationTrack", function(from, to) AnnotationTrack(range = from))

setAs("GRangesList", "AnnotationTrack", function(from, to) AnnotationTrack(range = from))

## Show ----------------------------------------------------------------------

## A helper function to plot general information about an AnnotationTrack
.annotationTrackInfo <- function(object) {
    msg <- sprintf(
        paste("| genome: %s\n| active chromosome: %s\n",
            "| annotation features: %s",
            sep = ""
        ),
        genome(object),
        chromosome(object),
        length(object)
    )
    addfeat <- length(object@range) - length(object)
    if (addfeat > 0) {
        msg <- c(msg, .addFeatInfo(object, addfeat), "Call chromosome(obj) <- 'chrId' to change the active chromosome")
    }
    return(paste(msg, collapse = "\n"))
}

## We have to show the name, genome and currently active chromosome, and, if more ranges are available on additional
## chromosomes some information about that

#' @describeIn AnnotationTrack-class Show method.
#' @export
setMethod("show", signature(object = "AnnotationTrack"), function(object) {
    cat(sprintf("AnnotationTrack '%s'\n%s\n", names(object), .annotationTrackInfo(object)))
})

#' @describeIn AnnotationTrack-class Show method.
#' @export
setMethod("show", signature(object = "ReferenceAnnotationTrack"), function(object) {
    .referenceTrackInfo(object, "ReferenceAnnotationTrack")
})