R/Gviz-methods.R

Defines functions sequenceTrackInfo referenceTrackInfo annotationTrackInfo addFeatInfo fillWithDefaults getBioColorIdeo roundedCap z2icol extreme buildArgsDf expLabel barsAndLabels boxes aggregator collapseAnnotation myFindOverlaps setAnn getAnn

##----------------------------------------------------------------------------------------------------------------------------
## Some rather general accessors to extract information from all kinds of GdObjects
##----------------------------------------------------------------------------------------------------------------------------
## Extract the full GRanges object from the range slot of an object inheriting from RangeTrack
setMethod("ranges", "RangeTrack", function(x) x@range)
setReplaceMethod("ranges", "RangeTrack", function(x, value) {
    x@range <- value
    return(x)})
setMethod("ranges", "GenomeAxisTrack", function(x) x@range)
setReplaceMethod("ranges", "GenomeAxisTrack", function(x, value) {
    x@range <- value
    return(x)})


## Extract the IRanges part of the GRanges object from the range slot of an object inheriting from RangeTrack
setMethod("range", "RangeTrack", function(x) ranges(x@range))
setMethod("range", "GenomeAxisTrack", function(x) ranges(x@range))

## seqnames, levels and infofrom the range track
setMethod("seqnames", "RangeTrack", function(x) as.character(seqnames(ranges(x))))
setMethod("seqnames", "SequenceDNAStringSetTrack", function(x) as.character(names(x@sequence)))
setMethod("seqnames", "SequenceBSgenomeTrack", function(x) as.character(seqnames(x@sequence)))
setMethod("seqlevels", "RangeTrack", function(x) unique(seqnames(x)))
setMethod("seqlevels", "SequenceDNAStringSetTrack", function(x) seqnames(x)[width(x@sequence)>0])
setMethod("seqlevels", "SequenceBSgenomeTrack", function(x) seqnames)
setMethod("seqinfo", "RangeTrack", function(x) table(seqnames(x)))

## Min and max ranges
setMethod("min", "RangeTrack", function(x) min(start(x)))
setMethod("max", "RangeTrack", function(x) max(end(x)))

## Extract start and end coordinates
setMethod("start", "RangeTrack", function(x) if(length(x)) as.integer(start(range(x))) else NULL)
setReplaceMethod("start", "RangeTrack", function(x, value) {
    start(x@range) <- value
    return(x)})
setReplaceMethod("start", "GenomeAxisTrack", function(x, value) {
    start(x@range) <- value
    return(x)})
setReplaceMethod("start", "IdeogramTrack", function(x, value) return(x))
setMethod("end", "RangeTrack", function(x) if(length(x)) as.integer(end(range(x))) else NULL)
setReplaceMethod("end", "RangeTrack", function(x, value) {
    end(x@range) <- value
    return(x)})
setReplaceMethod("end", "GenomeAxisTrack", function(x, value) {
    end(x@range) <- value
    return(x)})
setReplaceMethod("end", "IdeogramTrack", function(x, value) return(x))
setMethod("width", "RangeTrack", function(x) if(length(x)) as.integer(width(range(x))) else NULL)
setReplaceMethod("width", "RangeTrack", function(x, value) {
    width(x@range) <- value
    return(x)})
setReplaceMethod("width", "IdeogramTrack", function(x, value) return(x))
setMethod("start", "GenomeAxisTrack", function(x) if(length(x)) start(range(x)) else NULL)
setMethod("end", "GenomeAxisTrack", function(x) if(length(x)) end(range(x)) else NULL)
setMethod("width", "GenomeAxisTrack", function(x) if(length(x)) as.integer(width(range(x))) else NULL)
setMethod("start", "IdeogramTrack", function(x) NULL)
setMethod("start", "SequenceTrack", function(x) NULL)
setMethod("end", "IdeogramTrack", function(x) NULL)
setMethod("end", "SequenceTrack", function(x) NULL)
setMethod("width", "IdeogramTrack", function(x) NULL)
setMethod("width", "SequenceTrack", function(x) NULL)

## Return the number of individual annotation items (independent of any grouping) in a RangeTrack
setMethod("length", "RangeTrack", function(x) sum(seqnames(x) == chromosome(x)))
setMethod("length", "GenomeAxisTrack", function(x) length(ranges(x)))
setMethod("length", "IdeogramTrack", function(x) length(ranges(x)))
setMethod("length", "SequenceTrack", function(x)
          if(chromosome(x) %in% seqnames(x)) length(x@sequence[[chromosome(x)]]) else 0)
## setMethod("length", "ReferenceAnnotationTrack", function(x) 0)
## setMethod("length", "ReferenceGeneRegionTrack", function(x) 0)
## setMethod("length", "ReferenceDataTrack", function(x) 0)
setMethod("length", "HighlightTrack", function(x) length(x@trackList))
setMethod("length", "OverlayTrack", function(x) length(x@trackList))

## Extract the metadata columns from the GRanges object of an object inheriting from RangeTrack as a data.frame.
## For a DataTrack object these values are stored as a numeric matrix in the data slot, and we return this instead.
setMethod("values", "RangeTrack", function(x) as.data.frame(values(ranges(x))))
setMethod("values", "GenomeAxisTrack", function(x) as.data.frame(values(ranges(x))))
setMethod("values", "DataTrack", function(x, all=FALSE){
    if(sum(dim(x@data))==0) x@data else{
        sel <- if(all) rep(TRUE, ncol(x@data)) else seqnames(x) == chromosome(x)
        x@data[,sel, drop=FALSE]
    }
})
setMethod("values", "AlignmentsTrack", function(x) .dpOrDefault(x, ".__coverage"))
setReplaceMethod("values", "DataTrack", function(x, value){
    if(!is.matrix(value))
    {
        if(!is.numeric(value) || length(value) != length(x))
            stop("Invalid length of replacement vector.")
        if(!is.matrix(value) || !is.numeric(value) || ncol(value)!=length(x))
            stop("Dimensions of replacement value do not match.")
    }
    x@data <- value
    return(x)
})


## Extract a subsequence from a SequenceTrack. For performance reasons we restrict this to a maximum
## of one million nucleotides (which is already more than plenty...)
setMethod("subseq", "SequenceTrack", function(x, start=NA, end=NA, width=NA){
    padding <- "-"
    if(!is.na(start[1]+end[1]+width[1])){
        warning("All 'start', 'stop' and 'width' are provided, ignoring 'width'")
        width <- NA
    }
    ## We want start and end to be set if width is provided
    if(!is.na(width[1])){
        if(is.na(start) && is.na(end))
            stop("Two out of the three in 'start', 'end' and 'width' have to be provided")
        if(is.na(start))
            start <- end-width[1]+1
        if(is.na(end))
            end <- start+width[1]-1
    }
    w <- length(x)
    if(is.na(start))
        start <- 1
    if(w>0){
        if(is.na(end))
            end <- w
        rstart <- max(1, start[1], na.rm=TRUE)
        rend <- max(rstart, min(end[1], w, na.rm=TRUE))
    }else{
        if(is.na(end))
            end <- start
        rend <- end
        rstart <- start
    }
    if(rend<rstart)
        stop("'end' has to be bigger than 'start'")
    if((rend-rstart+1)>10e6)
        stop("Sequence is too big! Unable to extract")
    finalSeq <- rep(DNAString(padding), end-start+1)
    if(chromosome(x) %in% seqnames(x) && rend>rstart){
        chrSeq <- x@sequence[[chromosome(x)]]
        seq <- subseq(chrSeq, start=rstart, end=rend)
        if(is(x, "SequenceBSgenomeTrack")) seq <- unmasked(seq)
        subseq(finalSeq, ifelse(start<1, abs(start)+2, 1), width=rend-rstart+1) <- seq
    }
    if(is(x, "SequenceBSgenomeTrack") && chromosome(x) %in% seqnames(x))
        x@pointerCache[[chromosome(x)]] <- x@sequence[[chromosome(x)]]
    if(.dpOrDefault(x, "complement", FALSE))
        finalSeq <- complement(finalSeq)
    return(finalSeq)
})

setMethod("subseq", "ReferenceSequenceTrack", function(x, start=NA, end=NA, width=NA){
    ## We want start and end to be set if width is provided
    if(!is.na(width[1])){
        if(is.na(start) && is.na(end))
            stop("Two out of the three in 'start', 'end' and 'width' have to be provided")
        if(is.na(start))
            start <- end-width[1]+1
        if(is.na(end))
            end <- start+width[1]-1
    }
    x@sequence <- x@stream(file=x@reference, selection=GRanges(chromosome(x), ranges=IRanges(start, end)))
    return(callNextMethod())
})



## Set or extract the chromosome from a RangeTrack object
setMethod("chromosome", "GdObject", function(GdObject) return(NULL))
setMethod("chromosome", "RangeTrack", function(GdObject) GdObject@chromosome)
setMethod("chromosome", "SequenceTrack", function(GdObject) GdObject@chromosome)
setReplaceMethod("chromosome", "GdObject", function(GdObject, value){
    return(GdObject)
})
setReplaceMethod("chromosome", "RangeTrack", function(GdObject, value){
    GdObject@chromosome <- .chrName(value[1])
    return(GdObject)
})
setReplaceMethod("chromosome", "SequenceTrack", function(GdObject, value){
    GdObject@chromosome <- .chrName(value[1])
    return(GdObject)
})
setReplaceMethod("chromosome", "IdeogramTrack", function(GdObject, value){
    ## We have changed the class definition to include the bands for all chromosomes, but still want the old objects to work
    chromosome <- .chrName(value[1])
    if(.hasSlot(GdObject, "bandTable") && chromosome %in% as.character(GdObject@bandTable$chrom))
    {
        ranges <- GdObject@bandTable[GdObject@bandTable$chrom==chromosome,]
        bnames <- as.character(ranges$name)
        sel <- is.na(bnames)
        if(any(sel))
            bnames[sel] <- paste("band", seq_len(sum(sel)), sep="_")
        if(any(bnames == ""))
            bnames[bnames == ""] <- sprintf("band_%i", which(bnames == ""))
        ranges <- GRanges(seqnames=bnames, ranges=IRanges(start=ranges$chromStart, end=ranges$chromEnd),
                          name=bnames, type=ranges$gieStain)
        GdObject@range <- ranges
        GdObject@chromosome <- chromosome
        return(GdObject)
    }
    message("Updating chromosome band information")
    tmp <- IdeogramTrack(genome=genome(GdObject), chromosome=.chrName(value[1]), name=names(GdObject))
    displayPars(tmp) <- displayPars(GdObject, hideInternal=FALSE)
    return(tmp)
})

setReplaceMethod("chromosome", "AlignmentsTrack", function(GdObject, value){
    GdObject <- callNextMethod()
    if(!is.null(GdObject@referenceSequence))
        chromosome(GdObject@referenceSequence) <- value[1]
    return(GdObject)
})

setReplaceMethod("chromosome", "HighlightTrack", function(GdObject, value){
    GdObject@trackList <- lapply(GdObject@trackList, function(x){
        chromosome(x) <- value[1]
        x
    })
    GdObject@chromosome <- .chrName(value[1])
    return(GdObject)
})

setReplaceMethod("chromosome", "OverlayTrack", function(GdObject, value){
    GdObject@trackList <- lapply(GdObject@trackList, function(x){
        chromosome(x) <- value
        x
    })
    return(GdObject)
})


## Set or extract the genome from a RangeTrack object
setMethod("genome", "RangeTrack", function(x) x@genome)
setMethod("genome", "SequenceTrack", function(x) x@genome)
setReplaceMethod("genome", "GdObject", function(x, value){
    return(x)
})
setReplaceMethod("genome", "RangeTrack", function(x, value){
    x@genome <- value[1]
    genome(ranges(x)) <- as.vector(value[1])
    return(x)
})
setReplaceMethod("genome", "IdeogramTrack", function(x, value){
    if(genome(x)!=value)
        message("Updating chromosome band information")
    tmp <- IdeogramTrack(genome=value[1], chromosome=chromosome(x), name=names(x))
    displayPars(tmp) <- displayPars(x, hideInternal=FALSE)
    return(tmp)
})

## Set or extract the name slot of a GdObject
setMethod("names", "GdObject", function(x) x@name)
setReplaceMethod("names", signature("GdObject", "character"), function(x, value) {
    x@name <- value[1]
    return(x)
})

## Set or extract the strand information from a RangeTrack object
setMethod("strand", "RangeTrack", function(x) as.character(strand(ranges(x))))
setMethod("strand", "GenomeAxisTrack", function(x) as.character(strand(ranges(x))))
setMethod("strand", "DataTrack", function(x) x@strand)
setReplaceMethod("strand", "RangeTrack", function(x, value){
    if(length(value)!=1 && length(value)!=length(x))
        stop("Length of replacement value for the strand information does not match the ",
             "number of items in the track")
    r <- ranges(x)
    strand(r) <- value
    x@range <- r
    return(x)
})
setReplaceMethod("strand", "DataTrack", function(x, value){
    if(!is.character(value) && length(value)!=1 && !value %in% c("+", "-", "*"))
        stop("Invalid replacement value")
    x@strand <- value
    return(x)
 })

## Allow for subsetting of RangeTrack and DataTrack objects
setMethod("[", signature(x="RangeTrack"), function(x, i, j, ..., drop=TRUE) {
    x <- .deepCopyPars(x)
    x@range <- x@range[i,]
    return(x)})
setMethod("[", signature(x="StackedTrack"), function(x, i, j, ..., drop=TRUE) {
    x <- callNextMethod(x,i)
    x@stacks <- x@stacks[i]
    return(x)})
setMethod("[", signature(x="GenomeAxisTrack"), function(x, i, j, ..., drop=TRUE) {
    x <- .deepCopyPars(x)
    x@range <- x@range[i,]
    return(x)})
setMethod("[", signature(x="IdeogramTrack"), function(x, i, j, ..., drop=TRUE) return(x))
setMethod("[", signature(x="DataTrack"), function(x, i, j,  ..., drop=TRUE) {
    x <- .deepCopyPars(x)
    if(!missing(i))
    {
        x@data <- x@data[i,, drop=FALSE]
        displayPars(x) <- list(groups=as.vector(.dpOrDefault(x, "groups")[i]))
    }
    if(!missing(j))
    {
        x@range <- x@range[j,]
        if(ncol(x@data)>0)
            x@data <- x@data[,j, drop=FALSE]
    }
    return(x)})
setMethod("[", signature(x="AlignedReadTrack"), function(x, i, j, ..., drop=TRUE) {
    if(x@coverageOnly)
        stop("This AlignedReadTrack object contains coverage information only and can not be subset")
    x@range <- x@range[i,]
    x <- setCoverage(x)
    return(x)})

## Split a RangeTrack or DataTrack by a factor or character
setMethod("split", signature("RangeTrack"),
          definition=function(x, f, ...){
              rs <- split(ranges(x), factor(f))
              lapply(rs, function(y) {x@range <- y; return(x)})
          })
setMethod("split", signature("AlignedReadTrack"),
          definition=function(x, f, ...){
              if(x@coverageOnly)
                  stop("This AlignedReadTrack object contains coverage information only and can not be split")
              rs <- split(ranges(x), factor(f))
              lapply(rs, function(y) {x@range <- y; x <- setCoverage(x); return(x)})
		  })
setMethod("split", signature("DataTrack"),
		  definition=function(x, f, ...){
			  rs <- as.list(split(ranges(x), factor(f)))
			  ds <- split(t(values(x)), f)
			  nr <- nrow(values(x))
			  mapply(function(y, z) {x@range <- y; x@data <- matrix(z, nrow=nr, byrow=TRUE); return(x)}, rs, ds)
		  })

## Extract the coverage information
setMethod("coverage", signature("AlignedReadTrack"),
		definition=function(x, strand="*"){
			str <- c("+", "-", "*")[.strandName(strand, extended=TRUE)+1]
			return(if(!is.null(x@coverage[[str]])) x@coverage[[str]] else Rle())
		})

##----------------------------------------------------------------------------------------------------------------------------
## There are several levels of annotation information for most RangeTrack objects: individual features (e.g. exons, biotype),
## groups (e.g. transcripts) and even groups of groups (e.g. genes). Not all are relevant for all subclasses, however we want
## to have accessors and replacement methods for a clean interface.
##----------------------------------------------------------------------------------------------------------------------------
## Helper functions to extract or replace the various annotation data of a track.
##   o GdObject: the input GeneRegionTrack track object
##   o type: the annotation type, i.e., a metadata column in the GRanges object
##   o value: the replacement value, has to be of the same length as length(GdObject)
.getAnn <- function(GdObject, type) return(as.character(values(GdObject)[[type]]))
.setAnn <-  function(GdObject, value, type)
{
    v <- values(GdObject)
    if(length(value)>1 && length(value) != nrow(v))
        stop("The length of the replacement value for the '", type, "' annotation does not match the number ",
             "of features in the track.")
    v[[type]] <- value
    mcols(GdObject@range) <- v
    return(GdObject)
}

## Accessors to the gene-level annotation of a GeneRegionTrack. We actually need two methods here, one for the
## human-readable gene symbols and one for the actual gene id
setMethod("gene", signature(GdObject="GeneRegionTrack"), function(GdObject) .getAnn(GdObject, "gene"))
setMethod("symbol", signature(GdObject="GeneRegionTrack"), function(GdObject) .getAnn(GdObject, "symbol"))
setReplaceMethod("gene", signature("GeneRegionTrack", "character"), function(GdObject, value) .setAnn(GdObject, value, "gene"))
setReplaceMethod("symbol", signature("GeneRegionTrack", "character"), function(GdObject, value) .setAnn(GdObject, value, "symbol"))

## Accessors to the transcript-level annotation of a GeneRegionTrack.
setMethod("transcript", signature(GdObject="GeneRegionTrack"), function(GdObject) .getAnn(GdObject, "transcript"))
setReplaceMethod("transcript", signature("GeneRegionTrack", "character"),
                 function(GdObject, value) .setAnn(GdObject, value, "transcript"))

## Accessors to the exon-level annotation of a GeneRegionTrack
setMethod("exon", signature(GdObject="GeneRegionTrack"), function(GdObject) .getAnn(GdObject, "exon"))
setReplaceMethod("exon", signature("GeneRegionTrack", "character"), function(GdObject, value) .setAnn(GdObject, value, "exon"))

## Accessors to the biotype annotation of a RangeTrack.
setMethod("feature", signature(GdObject="RangeTrack"), function(GdObject) .getAnn(GdObject, "feature"))
setMethod("feature", signature(GdObject="DataTrack"), function(GdObject) NULL)
setReplaceMethod("feature", signature("RangeTrack", "character"), function(GdObject, value) .setAnn(GdObject, value, "feature"))
setReplaceMethod("feature", signature("DataTrack", "character"), function(GdObject, value) GdObject)

## Accessors to the grouping information of a AnnotationTrack and a GeneRegionTrack (for which this is essentially
## an alias to the transcript accessor.
setMethod("group", "AnnotationTrack", function(GdObject) .getAnn(GdObject, "group"))
setReplaceMethod("group", signature("AnnotationTrack", "character"), function(GdObject, value) .setAnn(GdObject, value, "group"))
setMethod("group", "GeneRegionTrack", function(GdObject) transcript(GdObject))
setReplaceMethod("group", signature("GeneRegionTrack", "character"),
                 function(GdObject, value) .setAnn(GdObject, value, "transcript"))
setMethod("group", "GdObject", function(GdObject) NULL)

## extract or replace the content of the imageMap slot
setMethod("imageMap", "GdObject", function(GdObject) GdObject@imageMap)
setReplaceMethod("imageMap", signature("GdObject", "ImageMapOrNULL"), function(GdObject, value){
    GdObject@imageMap <- value
    return(GdObject)})

## 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).
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)
})
setMethod("identifier", "GeneRegionTrack", function(GdObject, type=.dpOrDefault(GdObject, "transcriptAnnotation", "symbol")){
    if(is.logical(type)){
        type <- ifelse("symbol", "gene", type[1])
    }
    if(is.null(type)){
        type <- "symbol"
    }
    id <- switch(as.character(type),
                 "symbol"=symbol(GdObject),
                 "gene"=gene(GdObject),
                 "transcript"=transcript(GdObject),
                 "feature"=feature(GdObject),
                 "exon"=exon(GdObject),
                 "lowest"=exon(GdObject),
                 symbol(GdObject))
    id[is.na(id)] <- "NA"
    return(id)
})
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)
})
setReplaceMethod("identifier", c("GeneRegionTrack", "character"), function(GdObject, value){
    type <- .dpOrDefault(GdObject, "transcriptAnnotation", "symbol")
    switch(as.character(type),
           "symbol"=symbol(GdObject) <- value,
                 "gene"=gene(GdObject) <- value,
                 "transcript"=transcript(GdObject) <- value,
                 "feature"=feature(GdObject) <- value,
                 "exon"=exon(GdObject) <- value,
                 "lowest"=exon(GdObject) <- value,
                 symbol(GdObject) <- value)
    return(GdObject)
})
##----------------------------------------------------------------------------------------------------------------------------



##----------------------------------------------------------------------------------------------------------------------------
## Stacking controls what to do with overlapping annotation regions.
##----------------------------------------------------------------------------------------------------------------------------
setMethod("stacking", "StackedTrack", function(GdObject) GdObject@stacking)
setReplaceMethod("stacking", c("StackedTrack", "character"),
                 function(GdObject, value) {
                     pt <- getClass("StackedTrack")@prototype
                     if(!all(value %in% pt@stackingValues))
                         stop("Problem initializing StackedTrack,  need the following values for 'stacking':",
                              paste(pt@stackingValues, collapse=", "), "\n")
                     GdObject@stacking <- value
                     displayPars(GdObject) <- list(stacking=value)
                     return(GdObject)
                 })
##----------------------------------------------------------------------------------------------------------------------------



##----------------------------------------------------------------------------------------------------------------------------
## Recompute or return the stacking information for different types of StackedTrack objects. Stacking in needed when
## annotation regions in the objects overlap and when the stacking type is set to squish, full or pack. Since stacking
## can be dependent on the available space (if feature annotation is added) we need to be able to recompute this before
## we start the actual plotting. For the different sub-classes of StackedTracks we need different behaviour of setStacks:
##   o StackedTrack: there are no groups, so each feature can be treated separately
##   o AnnotationTrack and GeneRegionTrack: features can be grouped, in which case we have to avoid overlapping of the whole group region,
##      i.e, from the start of the first group item to the end of the last. In addition to grouping we have to factor in additional space
##      before each group if group/transcript annotation is enabled (gpar showId==TRUE). To do so we need to figure out the current
##      fontsize on the device, which means that a device already has to be open and the appropriate viewport has been
##      pushed to the stack. Hence we have to call setStacks immediately before the actual plotting.
## 'stacks' should return a vector of stacks, where each factor level of the vector indicates membeship
## to a particular stacking level. 'setStacks' returns the updated GdObject.
##----------------------------------------------------------------------------------------------------------------------------
setMethod("stacks", "StackedTrack",
          function(GdObject) if(length(GdObject@stacks)) GdObject@stacks else NULL)

setMethod("stacks", "AlignmentsTrack",
          function(GdObject) if(length(GdObject)) ranges(GdObject)$stack else 0)

setMethod("setStacks", "GdObject", function(GdObject, ...) GdObject)
setMethod("setStacks", "StackedTrack", function(GdObject, ...) {
    bins <- if(!.needsStacking(GdObject)) rep(1, length(GdObject)) else disjointBins(range(GdObject))
    GdObject@stacks <- bins
    return(GdObject)
})
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)
        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)
})
setMethod("setStacks", "HighlightTrack", function(GdObject, ...) {
     GdObject@trackList <- lapply(GdObject@trackList, setStacks, ...)
     return(GdObject)
 })
setMethod("setStacks", "OverlayTrack", function(GdObject, ...) {
     GdObject@trackList <- lapply(GdObject@trackList, setStacks, ...)
     return(GdObject)
 })

setMethod("setStacks", "AlignmentsTrack", function(GdObject, ...) {
    if(length(GdObject)){
        bins <- if(!.needsStacking(GdObject)) rep(1, length(GdObject)) else disjointBins(GdObject@stackRanges)
        GdObject@stacks <- bins
        ranges(GdObject)$stack <- bins[match(ranges(GdObject)$groupid, names(GdObject@stackRanges))]
    }
    return(GdObject)
})
##----------------------------------------------------------------------------------------------------------------------------



##----------------------------------------------------------------------------------------------------------------------------
## In some cases we don't need a range but rather a single position. Most of the time this is simply taking the geometric
## mean of the range. If numeric values are associated to positions we have to be able to extract those as well.
##----------------------------------------------------------------------------------------------------------------------------
## The geometric mean of annotation ranges
setMethod("position", signature("RangeTrack"), definition=function(GdObject, from=NULL, to=NULL, sort=FALSE, ...)
      {
          if(!is.null(from) && !is.null(to))
              GdObject <- subset(GdObject, from=from, to=to, sort=sort, ...)
          pos <- if(length(GdObject)) rowMeans(cbind(start(GdObject), end(GdObject))) else numeric()
          return(pos)
      })
setMethod("position", signature("IdeogramTrack"), definition=function(GdObject, ...) NULL)

## The numeric values of data tracks
setMethod("score", signature("DataTrack"), function(x, from=NULL, to=NULL, sort=FALSE, transformation=TRUE, ...)
      {
          if(!is.null(from) && !is.null(to))
              x <- subset(x, from=from, to=to, sort=sort, ...)
          vals <- values(x)
          ## apply data transformation if one is set up
          trans <- .dpOrDefault(x, "transformation")
          if(is.list(trans))
              trans <- trans[[1]]
          if(transformation && !is.null(trans))
          {
              if(!is.function(trans) || length(formals(trans))!=1L)
                  stop("gpar 'transformation' must be a function with a single argument")
              test <- trans(vals)
              if(!is.numeric(test) || !is.matrix(test) || !all(dim(test) == dim(vals)))
                  stop("The function in gpar 'transformation' results in invalid output.\n",
                       "It has to return a numeric matrix with the same dimensions as the input data.")
              vals <- test
          }
          return(vals)
      })



##----------------------------------------------------------------------------------------------------------------------------
## Before starting of the plotting operation there are a bunch of housekeeping task that should be performed on each
## track, and the mileage may vary between track types, hence we add a layer of abstraction here by using a method.
## Available arguments are:
##    o GdObject: the input track object
##    o chromosome: the currently active chromosome which may have to be set for a RangeTrack or a SequenceTrack object
##    o ...: additional arguments that are considered to be display parameters
##----------------------------------------------------------------------------------------------------------------------------
## For all track types we want to update the display parameters
setMethod("consolidateTrack", signature(GdObject="GdObject"), function(GdObject, alpha, ...) {
    pars <- list(...)
    pars <- pars[names(pars)!=""]
    pars[[".__hasAlphaSupport"]] <- alpha
    displayPars(GdObject) <- pars
    return(GdObject)
})
## For RangeTracks and SequenceTracks we want to set the chromosome
setMethod("consolidateTrack", signature(GdObject="RangeTrack"), function(GdObject, chromosome, ...) {
    if(!is.null(chromosome))
        chromosome(GdObject) <- chromosome
    GdObject <- callNextMethod(GdObject, ...)
    return(GdObject)
})
setMethod("consolidateTrack", signature(GdObject="SequenceTrack"), function(GdObject, chromosome, ...) {
    if(!is.null(chromosome))
        chromosome(GdObject) <- chromosome
    GdObject <- callNextMethod(GdObject, ...)
    return(GdObject)
})
## For StackedTracks we want to set the stacking (which could have been passed in as a display parameter)
setMethod("consolidateTrack", signature(GdObject="StackedTrack"), function(GdObject, ...) {
    GdObject <- callNextMethod()
    st <- .dpOrDefault(GdObject, "stacking")
    if(!is.null(st))
        stacking(GdObject) <- st
    return(GdObject)
})
## For AnnotationTracks we need to determine whether there is group label annotation or not, and we add this
## information as the internal display parameter '.__hasAnno'. We also precompute the grouped ranges together
## with optional labels in order to determine the correct plotting range later.
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)
})
## For a HighlightTrack we apply the method on each of the subtracks in the trackList slot
setMethod("consolidateTrack", signature(GdObject="HighlightTrack"), function(GdObject, chromosome, ...) {
    GdObject@trackList <- lapply(GdObject@trackList, consolidateTrack, chromosome=chromosome, ...)
    return(GdObject)
})
setMethod("consolidateTrack", signature(GdObject="OverlayTrack"), function(GdObject, chromosome, ...) {
    GdObject@trackList <- lapply(GdObject@trackList, consolidateTrack, chromosome=chromosome, ...)
    return(GdObject)
})
##----------------------------------------------------------------------------------------------------------------------------



##----------------------------------------------------------------------------------------------------------------------------
## There is a natural limit of what can be plotted as individual features caused by the maximum resolution of the device.
## Essentially no object can be smaller than the equivalent of a single pixel on the screen or whatever a pixel corresponds
## to on other devices.
## Thus we need to adapt the all objects in the track to the current resolution. Things that are closer than a certain limit
## distance can be collapsed (if that is desired), things that are smaller than 1 pixel can be blown up to a minimum size.
## All collapseTrack methods should always return a GdObject instance of similar type as the input to allow us to keep the
## collpasing step optional and for all the downstream plotting operations to still work. The internal (mostly the GRanges)
## objects can be modified to whatever is desired. Please note that this method is called after the drawing viewport has been
## pushed, and hence all the coordinate systems are already in place.
## Available arguments are:
##    o GdObject: the input AnnotationTrack or GeneRegionTrack track object
##    o min.width: the minimum width in pixel, everything that's smaller will be expanded to this size.
##    o min.distance: the minimum distance between two features in pixels below which to start collapsing track items.
##    o collapse: logical, collapse overlapping items into a single meta-item.
##    o diff: the equivalent of 1 pixel in the native coordinate system.
##    o xrange: the data range on the x axis. Can be used for some preliminary subsetting to speed things up
##----------------------------------------------------------------------------------------------------------------------------
## 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
        newVals <- rbind(newVals, as.data.frame(t(sapply(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"]]))
            }
        })), 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)
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(sapply(xs, function(y) abs(diff(as.numeric(as.data.frame(ranges(range(y)))[,c("start", "end")])))))]]
                             }))
                             GdObject@range <- r
                         },
                         "shortest"={
                             r <- unlist(endoapply(split(r, gene(GdObject)), function(x){
                                 xs <- split(x, x$transcript)
                                 xs[[which.min(sapply(xs, function(y) abs(diff(as.numeric(as.data.frame(ranges(range(y)))[,c("start", "end")])))))]]
                             }))
                             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)})



.aggregator <- function(GdObject)
{
    agFun <- .dpOrDefault(GdObject, "aggregation", "mean")
    if(is.list(agFun))
        agFun <- agFun[[1]]
    fun <- if(is.character(agFun))
    {
        switch(agFun,
               "mean"=rowMeans,
               "sum"=rowSums,
               "median"=rowMedians,
               "extreme"=function(x) apply(x, 1, .extreme),
               "min"=rowMin,
               "max"=rowMax,
               rowMeans)
    } else {
        if(is.function(agFun)){
            function(x) apply(x, 1, agFun)
        }else stop("display parameter 'aggregation' has to be a function or a character",
                   "scalar in c('mean', 'median', 'sum', 'extreme')")
    }
    return(fun)
}



## For DataTracks we want to collapse data values using the aggregation function provided by calling .aggregator().
## In addition values can be aggregated over fixed window slices when gpar 'window' is not NULL, and using a sliding
## window approach when 'window' == -1
setMethod("collapseTrack", signature(GdObject="DataTrack"), function(GdObject, diff=.pxResolution(coord="x"), xrange) {
    if(!length(GdObject))
        return(GdObject)
    ## first the data transformation if needed
    values(GdObject) <- score(GdObject)
    collapse <- .dpOrDefault(GdObject, "collapse", FALSE)
    min.width <- .dpOrDefault(GdObject, "min.width", 2)
    min.distance <- max(0, .dpOrDefault(GdObject, "min.distance", 0))
    ## When an averaging window has been set, split the data up into these average chunks
    window <- .dpOrDefault(GdObject, "window")
    windowSize <- .dpOrDefault(GdObject, "windowSize")
    if(!is.null(window) || collapse)
        GdObject <- GdObject[,order(range(GdObject))]
    r <- ranges(GdObject)
    drange <- c(floor(xrange[1]), ceiling(xrange[2]))
    if(!is.null(window))
    {
        rr <-  if(is(r, "GRanges")) ranges(r) else r
        fw <- FALSE
        if(window=="auto")
            window <- min(ncol(values(GdObject)), 1000,
                          ceiling(width(range(rr))/(min.width*diff)))
        if(window=="fixed"){
            fw <- TRUE
            window <- 100
        }
        if(!is.numeric(window) || length(window)!=1L)
            stop("gpar 'window' must be a numeric scalar")
        window <- as.integer(window)
        sc <- values(GdObject)
        agFun <- .aggregator(GdObject)
        if(window==1)
        {
            sc <- matrix(agFun(sc), ncol=1)
            rtmp  <- IRanges(start=max(1,drange[1]), end=max(1, drange[2]-1))
            r <- if(is(r, "GRanges"))  GRanges(seqnames=seqnames(r)[1], ranges=rtmp) else rtmp
        } else if(window<1){
            if(is.null(windowSize))
                windowSize <- (max(GdObject)-min(GdObject))/100
            if(windowSize %% 2 !=1)
                windowSize <- windowSize+1
            rm <- vector("integer", width(range(range(GdObject))))
            ind <- unlist(mapply(function(x, y) x:y, start(GdObject), end(GdObject)))-min(GdObject)+1
            rm[ind] <- rep(sc[1,], width(GdObject))
            runwin <- suppressWarnings(runmean(Rle(as.numeric(rm)), k=windowSize, endrule="constant"))
            seqSel <- findRun(as.integer(position(GdObject))-min(GdObject)+1, runwin)
            newDat <- matrix(runValue(runwin)[seqSel], nrow=1)
            if(nrow(sc)>1)
            {
                newDat <- rbind(newDat,
                                matrix(sapply(2:nrow(sc), function(x) {
                                    rm[ind] <- rep(sc[x,], width(GdObject))
                                    suppressWarnings(runValue(runmean(Rle(as.numeric(rm)), k=windowSize, endrule="constant", na.rm=TRUE)))[seqSel]}),
                                       nrow=nrow(sc)-1, byrow=TRUE))
            }
            sc <- newDat
        } else {
            if(!is.null(window) && window > diff(drange))
                window <- diff(drange)
            if(!fw || is.null(windowSize)){
                windowSize <-  diff(drange) %/% window
            }else{
                window <- max(1, diff(drange) %/% windowSize)
            }
            remain <-  (diff(drange) - (window * windowSize))/2
            ir <- IRanges(start=seq(from=drange[1]+remain, to=drange[2]-remain-windowSize, length.out=window),
                          width=windowSize)
            if(remain>0)
                ir <- c(IRanges(start=drange[1], width=ceiling(remain)), ir,
                        IRanges(start=drange[2]-ceiling(remain), width=ceiling(remain)))
            ol <- as.matrix(findOverlaps(ir, rr))
            scn <- sapply(split(ol[,2], ol[,1]), function(i) agFun(sc[,i,drop=FALSE]), USE.NAMES=FALSE)
            if(is.null(dim(scn)))
                scn <- matrix(scn, nrow=nrow(values(GdObject)), dimnames=list(NULL, as.character(unique(ol[,1]))))
            sc <- matrix(NA, ncol=length(ir), nrow=nrow(scn))
            sc[, as.integer(colnames(scn))] <- scn
            r <- if(is(r, "GRanges")) GRanges(seqnames=chromosome(GdObject), ranges=ir,
                                              strand=unique(as.character(strand(GdObject)))) else ir
        }
        GdObject@range <- r
        GdObject@data <- sc
    }
    ## If groups need to be averaged we have to do it here
    groups <- .dpOrDefault(GdObject, "groups")
    if(!is.null(groups) && .dpOrDefault(GdObject, "aggregateGroups", FALSE)){
        if(!is.factor(groups))
            groups <- factor(groups)
        agFun <- .aggregator(GdObject)
        dat <- values(GdObject)
        rownames(dat) <- groups
        datNew <- matrix(t(sapply(levels(groups), function(x) agFun(t(dat[groups==x,,drop=FALSE])), USE.NAMES=FALSE)), nrow=nlevels(groups))
        GdObject@data <- datNew
        displayPars(GdObject) <- list(groups=levels(groups))
    }
    ## Compute native coordinate equivalent to 1 pixel and resize
    r <- .resize(r, min.width, diff)
    ## Collapse overlapping ranges (less than minXDist space between them) including the associated attributes using
    ## "|" as separator. For both "strand" and "feature" we take the first available entry, which is not optimal but
    ## seems to be the sanest thing to do here...
    if(collapse)
    {
        minXDist <- min.distance*diff
        rr <- if(is(r, "GRanges")) ranges(r) else r
        if(minXDist<1)
        {
            ## We have to fake smaller ranges because reduce will merge also neigbouring ranges
            width(rr) <- width(rr)-1
            rr <- reduce(rr, min.gapwidth=minXDist)
            width(rr) <- width(rr)+1
        } else {
            rr <- reduce(r, min.gapwidth=minXDist)
        }
        sc <- values(GdObject)
        if(length(rr)==1){
            r <- GRanges(seqnames=1, strand=strand(GdObject)[1], ranges=rr)
            GdObject@range <- r
            GdObject@data <- matrix(rowMeans(sc, na.rm=TRUE), ncol=1)
        } else if(length(rr) < length(r)){
            startInd <- sort(unique(sapply(start(rr), function(x) which(start(r)==x))))
            st <- strand(GdObject)
            startInd <- if(tail(startInd,1) == length(r)) c(startInd, length(r)+1) else c(startInd, length(r))
            vsplit <- split(t(as.data.frame(sc, stringsAsFactors=FALSE)), cut(seq_len(length(r)), startInd, iclude.lowest=TRUE, right=FALSE))
            agFun <- .dpOrDefault(GdObject, "aggregation", "mean")
            if(is.list(agFun))
                agFun <- agFun[[1]]
            newScore <- if(is.character(agFun)){
                switch(agFun, "mean"=sapply(vsplit, function(x) rowMeans(matrix(x, nrow=nrow(sc), byrow=TRUE), na.rm=TRUE), USE.NAMES=FALSE),
                       "sum"=sapply(vsplit, function(x) rowSums(matrix(x, nrow=nrow(sc), byrow=TRUE), na.rm=TRUE), USE.NAMES=FALSE),
                       "median"=sapply(vsplit, function(x) rowMedians(matrix(x, nrow=nrow(sc), byrow=TRUE), na.rm=TRUE), USE.NAMES=FALSE),
                       sapply(vsplit, function(x) rowMeans(matrix(x, nrow=nrow(sc), byrow=TRUE), na.rm=TRUE), USE.NAMES=FALSE))
            } else {
                if(is.function(agFun)){
                    sapply(vsplit, function(x) apply(matrix(x, nrow=nrow(sc), byrow=TRUE), 1, function(y) agFun(y)[1]), USE.NAMES=FALSE)
                } else stop("display parameter 'aggregation' has to be a function or a character ", "scalar in c('mean', 'median', 'sum')")
            }
            r <- GRanges(seqnames=seq_len(length(rr)), strand=st, ranges=rr)
            GdObject@data <- newScore
            GdObject@range <- r
        }
    }
    ## Reconstruct the RangedData object and return
    GdObject@range <- r
    return(GdObject)
})



## For a GenomeAxisTrack all we need to do is collapse the optional ranges
setMethod("collapseTrack", signature(GdObject="GenomeAxisTrack"),
          function(GdObject, min.width=1, min.distance=0, collapse=TRUE, diff=.pxResolution(coord="x"), xrange) {

              ## Collapse overlapping ranges (less than minXDist space between them) including the associated attributes using
              ## "|" as separator. For both "strand" and "feature" we take the first available entry, which is not optimal but
              ## seems to be the sanest thing to do here...
              if(collapse)
              {
                  GdObject <- GdObject[order(range(GdObject))]
                  r <- ranges(GdObject)
                  minXDist <- min.distance*diff
                  r <- reduce(r, min.gapwidth=minXDist)
              }
              r <- .resize(r, min.width, diff)
              GdObject@range <- r
              return(GdObject)})
##----------------------------------------------------------------------------------------------------------------------------



##----------------------------------------------------------------------------------------------------------------------------
## Truncate a GdObject and sort by coordinates if necessary.
##----------------------------------------------------------------------------------------------------------------------------
## The default is not to clip at all
setMethod("subset", signature(x="GdObject"), function(x, ...) x)

## For normal ranges we clip everything outside of the boundaries (keeping one extra item left and right
## in order to assure continuation)
setMethod("subset", signature(x="RangeTrack"), function(x, from=NULL, to=NULL, sort=FALSE, drop=TRUE, use.defaults=TRUE, ...){
    ## Not needed anymore...
    ## Subset to a single chromosome first
    if(drop){
        csel <- seqnames(x) != chromosome(x)
        if(any(csel))
            x <- x[,!csel]
    }
    if(!length(x))
        return(x)
    ranges <- if(use.defaults) .defaultRange(x, from=from, to=to) else c(from=ifelse(is.null(from), -Inf, from), to=ifelse(is.null(to), Inf, to))
    lsel <- end(x) < ranges["from"]
    if(any(lsel))
        lsel[max(0, max(which(lsel))-1)] <- FALSE
    rsel <- start(x) > ranges["to"]
    if(any(rsel))
        rsel[min(length(x), min(which(rsel))+1)] <- FALSE
    if(any(lsel) || any(rsel))
        x <- x[!(lsel | rsel),]
    if(sort)
        x <- x[order(range(x)),]
    return(x)
})


## For highlight track we just apply the method for all the subtracks in the tracklList and finally the RangeTrack method
setMethod("subset", signature(x="HighlightTrack"), function(x, ...){
    x@trackList <- lapply(x@trackList, subset, ...)
    x <- callNextMethod()
    return(x)
})

setMethod("subset", signature(x="OverlayTrack"), function(x, ...){
    x@trackList <- lapply(x@trackList, subset, ...)
    return(x)
})

## For DataTracks we cut exactly, and also reduce to the current chromosome unless told explicitely not to
setMethod("subset", signature(x="DataTrack"), function(x, from=NULL, to=NULL, sort=FALSE, drop=TRUE, use.defaults=TRUE, ...){
    ## Subset to a single chromosome first
    if(drop){
        csel <- seqnames(x) != chromosome(x)
        if(any(csel))
            x <- x[,!csel]
    }
    if(!length(x))
        return(x)
    ranges <- if(use.defaults) .defaultRange(x, from=from, to=to) else c(from=ifelse(is.null(from), -Inf, from), to=ifelse(is.null(to), Inf, to))
    x <- x[,start(x)>=ranges["from"] & end(x)<=ranges["to"]]
    if(sort)
        x <- x[,order(range(x))]
    return(x)
})

## ReferenceDataTracks need to stream the data from file and then pass the results on to the next method
setMethod("subset", signature(x="ReferenceDataTrack"), 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 ReferenceDataTrack")
    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))){
        vals <- x@stream(x@reference, subRegion)
        x@range <- vals
        mcols(x@range) <- NULL
        x@data <- .prepareDtData(if(ncol(values(vals))) as.data.frame(values(vals)) else matrix(nrow=0, ncol=0), length(vals))
        chromosome(x) <- chromosome[1]
    }
    return(callNextMethod(x=x, from=from, to=to, drop=FALSE, ...))
})

## Only recompute the stacks here
setMethod("subset", signature(x="StackedTrack"), function(x, from=NULL, to=NULL, sort=FALSE, stacks=FALSE, ...){
    x <- callNextMethod(x=x, from=from, to=to, sort=sort)
    if(stacks)
        x <- setStacks(x)
    return(x)
})

## 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.
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
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, ...))
})

## FIXME: Still needs to be implemented
setMethod("subset", signature(x="ReferenceGeneRegionTrack"), function(x, ...){
    warning("ReferenceGeneRegionTrack objects are not supported yet.")
    return(callNextMethod())
})

setMethod("subset", signature(x="BiomartGeneRegionTrack"), function(x, from, to, chromosome, use.defaults=TRUE, ...){
   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(ranges["from"] < x@start || ranges["to"] > x@end){
       x@start <- ranges["from"] - 10000
       x@end <- ranges["to"] + 10000
       ranges(x) <- .cacheMartData(x, .chrName(chromosome))
   }
   return(callNextMethod(x=x, from=ranges["from"], to=ranges["to"], use.defaults=FALSE, ...))
})


## For the axis track we may have to clip the highlight ranges on the axis.
setMethod("subset", signature(x="GenomeAxisTrack"), function(x, from=NULL, to=NULL, sort=FALSE, ...){
    if(!length(x))
        return(x)
    ranges <- .defaultRange(x, from=from, to=to)
    lsel <- end(x) < ranges["from"]
    rsel <- start(x) > ranges["to"]
    x <- x[!(lsel | rsel),]
    if(sort)
        x <- x[order(range(x)),]
    return(x)
})

## If the object only stores coverage we subset that, otherwise we can use the RangeTrack method
setMethod("subset", signature(x="AlignedReadTrack"), function(x, from=NULL, to=NULL, sort=FALSE, stacks=FALSE, ...){
    if(x@coverageOnly)
    {
        if(is.null(from))
            from <- min(unlist(lapply(x@coverage, function(y) if(length(y)) min(start(y)))))
        if(is.null(to))
            to <- max(unlist(lapply(x@coverage, function(y) if(length(y)) max(start(y)))))
        x@coverage <- lapply(x@coverage, function(y){runValue(y)[end(y)<from | start(y)>to] <- 0; y})
        x@coverage <- lapply(x@coverage, function(y){ if (length(y) < to) y <- c(y, Rle(0, to-length(y))); y})
        ##
        ##from <- min(unlist(lapply([email protected], function(y) if (length(y)) head(start(y), 2)[2])))
        if (max(unlist(lapply(x@coverage, function(y) {length(runLength(y)[runValue(y)!=0])}))))
        {
        from <- min(unlist(lapply(x@coverage, function(y) if(length(y)) head(start(y)[runValue(y)!=0],1))))
        to <- max(unlist(lapply(x@coverage, function(y) if(length(y)) tail(end(y),2)[1])))
        }
        x@range <- GRanges(ranges=IRanges(start=from, end=to), strand=names(x@coverage), seqnames=x@chromosome)
    }else{
        x <- callNextMethod(x=x, from=from, to=to, sort=sort, stacks=stacks)
    }
    return(x)
})

## AlignmentTracks can be subset by using the information in the stackRanges slot, but for the actual reads we need to make sure that
## we keep all the bits that belong to a given group. We still want to record the requested ranges in the internal '.__plottingRange'
## display parameter.
setMethod("subset", signature(x="AlignmentsTrack"), function(x, from=NULL, to=NULL, 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, otherwise we subset the stackRanges and keep all group items
        ranges <- if(use.defaults) .defaultRange(x, from=from, to=to) else c(from=ifelse(is.null(from), min(start(x@stackRanges))-1, from),
                                                                             to=ifelse(is.null(to), max(end(x@stackRanges))+1, to))
        displayPars(x) <- list(".__plottingRange"=ranges)
        sr <- subsetByOverlaps(x@stackRanges, GRanges(seqnames=chromosome(x)[1], ranges=IRanges(start=ranges["from"], end=ranges["to"])))
        if(length(sr) < length(x@stackRanges)){
            x@stackRanges <- sr
            x@range <- x@range[x@range$groupid %in% names(sr)]
            if(stacks)
                x <- setStacks(x)
        }
    }
    return(x)
})

## ReferenceAlignmentsTracks need to stream the data from file and then pass the results on to the next method
setMethod("subset", signature(x="ReferenceAlignmentsTrack"), 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)
    }else{
        chromosome <- .chrName(chromosome)
    }
    subRegion <- GRanges(seqnames=chromosome[1], ranges=IRanges(start=from, end=to))
    oldRange <- .dpOrDefault(x, ".__plottingRange")
    isIn <- length(ranges(x)) != 0 && !is.null(oldRange) && ranges(subRegion) %within% IRanges(oldRange["from"], oldRange["to"])
    if(!isIn){
        cMap <- .resolveColMapping(x@stream(x@reference, subRegion), x@args, x@mapping)
        seqs <- cMap$data$seq
        cMap$data$seq <- NULL
        range <- .computeAlignments(.buildRange(cMap$data, args=cMap$args, defaults=x@defaults, trackType="AnnotationTrack"))
        ranges(x) <- range$range
        x@stackRanges <- range$stackRanges
        x@stacks <- range$stacks
        x@sequences <- seqs
    }else{
       x@sequences <- subseq(x@sequences, start=from-(oldRange["from"]-1), width=to-from+1)
    }
    chromosome(x) <- chromosome[1]
    return(callNextMethod(x=x, from=from, to=to, drop=FALSE, ...))
})
##----------------------------------------------------------------------------------------------------------------------------




##----------------------------------------------------------------------------------------------------------------------------
## The indivdual bits and pieces of a Gviz plot are all drawn by separate renderers. Currently, those are a y-axis,
## a grid, and the actual track panel.
##----------------------------------------------------------------------------------------------------------------------------
## For certain GdObject subclasses we may want to draw a y-axis. For others an axis is meaningless, and the default function
## will return NULL without plotting anything.
setMethod("drawAxis", signature(GdObject="GdObject"), function(GdObject, ...) return(NULL))

setMethod("drawAxis", signature(GdObject="DataTrack"), function(GdObject, ...) {
    if(as.logical(.dpOrDefault(GdObject, "legend", FALSE)) && !is.null(.dpOrDefault(GdObject, ".__groupLevels"))){
        pushViewport(viewport(y=1, height=unit(1, "npc") - unit(.dpOrDefault(GdObject, ".__verticalSpace"), "inches"),
                              just=c(0.5, 1)))
        on.exit(popViewport(1))
    }
    type <- match.arg(.dpOrDefault(GdObject, "type", "p"), .PLOT_TYPES, several.ok=TRUE)
    isOnlyHoriz <- length(setdiff(type, "horizon")) == 0
    if(!isOnlyHoriz && .dpOrDefault(GdObject, "showAxis", TRUE)) {
	callNextMethod()
    } else {
        if(.dpOrDefault(GdObject, "showSampleNames", FALSE)){
            groups <- .dpOrDefault(GdObject, "groups")
            sn <- if(is.null(groups)) rownames(values(GdObject)) else rev(unlist(split(rownames(values(GdObject)), factor(groups))))
            cex.sn <- .dpOrDefault(GdObject, "cex.sampleNames", .dpOrDefault(GdObject, "cex.axis", 1))
            col.cn <- .dpOrDefault(GdObject, "col.sampleNames", "white")
            wd <- max(as.numeric(convertWidth(stringWidth(sn) + unit(10, "points"), "npc"))) * cex.sn
            samNames <- viewport(x=1, width=wd, just=1, yscale=c(-0.05, 1.05))
            pushViewport(samNames)
            nr <- nrow(values(GdObject))
            if(nr > 1){
                yy <- head(seq(0.05, 0.95, len=nr+1), -1)
                yy <- yy + diff(yy)[[1]]/2
            }else{
                yy <- 0.5
            }
            grid.text(x=rep(0.5, nr), y=yy, label=rev(sn), just=0.5, gp=gpar(cex=cex.sn, col=col.cn))
            popViewport(1)
        }
    }
})

setMethod("drawAxis", signature(GdObject="NumericTrack"), function(GdObject, from, to, ...) {
    type <- match.arg(.dpOrDefault(GdObject, "type", "p"), .PLOT_TYPES, several.ok=TRUE)
    yvals <- values(GdObject)
    ylim <- .dpOrDefault(GdObject, "ylim", if(!is.null(yvals) && length(yvals))
                         range(yvals, na.rm=TRUE, finite=TRUE) else c(-1,1))
    if(diff(ylim)==0)
        ylim <- ylim+c(-1,1)
    hSpaceAvail <- vpLocation()$isize["width"]/6
    yscale <- extendrange(r=ylim, f=0.05)
    col <- .dpOrDefault(GdObject, "col.axis", "white")
    acex <- .dpOrDefault(GdObject, "cex.axis")
    acol <- .dpOrDefault(GdObject, "col.axis", "white")
    at <- pretty(yscale)
    at <- at[at>=sort(ylim)[1] & at<=sort(ylim)[2]]
    if(is.null(acex))
    {
        vSpaceNeeded <- max(as.numeric(convertWidth(stringHeight(at), "inches")))*length(at)*1.5
        hSpaceNeeded <- max(as.numeric(convertWidth(stringWidth(at), "inches")))
        vSpaceAvail <- abs(diff(range(at)))/abs(diff(yscale))*vpLocation()$isize["height"]
        acex <- max(0.6, min(vSpaceAvail/vSpaceNeeded, hSpaceAvail/hSpaceNeeded))
    }
    nlevs <- max(1, nlevels(factor(.dpOrDefault(GdObject, "groups"))))
    if(type %in% c("heatmap", "horizon") && .dpOrDefault(GdObject, "showSampleNames", FALSE)){
        groups <- .dpOrDefault(GdObject, "groups")
        sn <- if(is.null(groups)) rownames(values(GdObject)) else rev(unlist(split(rownames(values(GdObject)), factor(groups))))
        cex.sn <- .dpOrDefault(GdObject, "cex.sampleNames", acex)
        col.cn <- .dpOrDefault(GdObject, "col.sampleNames", "white")
        wd <- max(as.numeric(convertWidth(stringWidth(sn) + unit(10, "points"), "npc"))) * cex.sn
        samNames <- viewport(x=1, width=wd, just=1, yscale=c(-0.05, 1.05))
        pushViewport(samNames)
        nr <- nrow(values(GdObject))
        if(nr > 1){
            yy <- head(seq(0.05, 0.95, len=nr+1), -1)
            yy <- yy + diff(yy)[[1]]/2
        }else{
            yy <- 0.5
        }
        grid.text(x=rep(0.5, nr), y=yy, label=rev(sn), just=0.5, gp=gpar(cex=cex.sn, col=col.cn))
        popViewport(1)
        samAxis <- viewport(x=1-wd, width=1-wd, just=1)
        pushViewport(samAxis)
        on.exit(popViewport(1))
    }
    ## if any of the types are gradient or heatmap we want the gradient scale
    if(any(type %in% c("gradient", "heatmap")) && .dpOrDefault(GdObject, "showColorBar", TRUE)){
        ## viewport to hold the color strip
        shift <- ifelse(all(type %in% c("gradient", "heatmap")), 1, 0)
        pcols <- .getPlottingFeatures(GdObject)
        ncolor <- .dpOrDefault(GdObject, "ncolor", 100)
        vpAxisCont <- viewport(x=unit(1, "npc")-unit(2-shift, "points"), width=unit(1, "npc")-unit(2-shift, "points"), just=1)
        pushViewport(vpAxisCont)
        for(i in seq_len(nlevs)){
            ## create color palette
	    cr <- c("white", pcols$col[i])
	    if(nlevs<2)
		cr <- .dpOrDefault(GdObject, "gradient", cr)
            palette <- colorRampPalette(cr)(ncolor+5)[-(1:5)]
            pshift <- ifelse(i==nlevs, 1-shift, 0)
            vpTitleAxis <- viewport(x=unit(1, "npc")-unit(4*(i-1), "points"), width=unit(4+pshift, "points"),
                                    yscale=yscale, just=1)
            pushViewport(vpTitleAxis)
            ## draw a rectangle for each color
            if(all(type %in% c("gradient", "heatmap"))){
                if(i==nlevs)
                    suppressWarnings(grid.yaxis(gp=gpar(col=acol, cex=acex), at=at))
                grid.rect(y=unit(seq(ylim[1],ylim[2],length.out=ncolor+1),"native")[-(ncolor+1)], x=unit(0, "npc")-unit(1, "points"),
                          width=1, height=1/ncolor, gp=gpar(fill=palette, lty=0), just=c("left","bottom"))
            }else{
                grid.rect(y=unit(seq(ylim[1],ylim[2],length.out=ncolor+1),"native")[-(ncolor+1)], x=0,
                          width=1, height=1/ncolor, gp=gpar(fill=palette, lty=0), just=c("left","bottom"))
                if(i==nlevs){
                    suppressWarnings(grid.yaxis(gp=gpar(col=acol, cex=acex), at=at))
                    grid.lines(x=c(0,0), y=ylim, gp=gpar(col=acol), default.units="native")
                }
            }
            popViewport(1)
        }
	popViewport(1)
    } else {
	vpTitleAxis <- viewport(x=0.95, width=0.2, yscale=yscale, just=0)
	pushViewport(vpTitleAxis)
	suppressWarnings(grid.yaxis(gp=gpar(col=acol, cex=acex), at=at))
	grid.lines(x=c(0,0), y=ylim, gp=gpar(col=acol), default.units="native")
	popViewport(1)
    }
})



setMethod("drawAxis", signature(GdObject="AlignmentsTrack"), function(GdObject, ...) {
    type <- match.arg(.dpOrDefault(GdObject, "type", .ALIGNMENT_TYPES), .ALIGNMENT_TYPES, several.ok=TRUE)
    if("coverage" %in% type){
        yvals <- values(GdObject)
        ylim <- .dpOrDefault(GdObject, "ylim", if(!is.null(yvals) && length(yvals))
                             range(yvals, na.rm=TRUE, finite=TRUE) else c(-1,1))
        if(diff(ylim)==0)
            ylim <- ylim+c(-1,1)
        hSpaceAvail <- vpLocation()$isize["width"]/6
        yscale <- c(if(is.null(.dpOrDefault(GdObject, "transformation"))) 0 else min(ylim), max(ylim) + diff(range(ylim)) * 0.05)
        col <- .dpOrDefault(GdObject, "col.axis", "white")
        acex <- .dpOrDefault(GdObject, "cex.axis")
        acol <- .dpOrDefault(GdObject, "col.axis", "white")
        at <- pretty(yscale)
        at <- at[at>=sort(ylim)[1] & at<=sort(ylim)[2]]
        covHeight <- .dpOrDefault(GdObject, ".__coverageHeight", c(npc=0, points=0))
        covSpace <- .dpOrDefault(GdObject, ".__coverageSpace", 0)
        pushViewport(viewport(y=1-(covHeight["npc"] + covSpace), height=covHeight["npc"], just=c(0.5, 0)))
        if(is.null(acex))
        {
            vSpaceNeeded <- max(as.numeric(convertWidth(stringHeight(at), "inches")))*length(at)*1.5
            hSpaceNeeded <- max(as.numeric(convertWidth(stringWidth(at), "inches")))
            vSpaceAvail <- abs(diff(range(at)))/abs(diff(yscale))*vpLocation()$isize["height"]
            acex <- max(0.6, min(vSpaceAvail/vSpaceNeeded, hSpaceAvail/hSpaceNeeded))
        }
        vpTitleAxis <- viewport(x=0.95, width=0.2, yscale=yscale, just=0)
        pushViewport(vpTitleAxis)
        suppressWarnings(grid.yaxis(gp=gpar(col=acol, cex=acex), at=at))
        grid.lines(x=c(0,0), y=ylim, gp=gpar(col=acol), default.units="native")
        popViewport(2)
    } else {
        covHeight <- c(npc=0, points=0)
        covSpace <- 0
    }
    if("sashimi" %in% type){
        sash <- .dpOrDefault(GdObject, ".__sashimi", list(x = numeric(), y = numeric(), id = integer(), score = numeric()))
        yscale <- if (length(sash$y)) c(-(max(sash$y) + diff(range(sash$y)) * 0.05), 0) else c(-1, 0)
        ylim <- if (length(sash$y)) c(-max(sash$y), yscale[1]+max(sash$y)) else c(-1,0)
        hSpaceAvail <- vpLocation()$isize["width"]/6
        col <- .dpOrDefault(GdObject, "col.axis", "white")
        acex <- .dpOrDefault(GdObject, "cex.axis")
        acol <- .dpOrDefault(GdObject, "col.axis", "white")
        labs <- if (length(sash$score)) pretty(c(1, sash$score)) else pretty(c(1, .dpOrDefault(GdObject, ".__sashimiScore", 10)))
        at <- seq(ylim[1], ylim[2], length.out=length(labs))
        sashHeight <- .dpOrDefault(GdObject, ".__sashimiHeight", c(npc=0, points=0))
        sashSpace <- .dpOrDefault(GdObject, ".__sashimiSpace", 0)
        pushViewport(viewport(y=1-(sashHeight["npc"] + sashSpace + covHeight["npc"] + covSpace),
                              height=sashHeight["npc"], just=c(0.5, 0)))
        if(is.null(acex))
        {
            vSpaceNeeded <- max(as.numeric(convertWidth(stringHeight(labs), "inches")))*length(at)*1.5
            hSpaceNeeded <- max(as.numeric(convertWidth(stringWidth(labs), "inches")))
            vSpaceAvail <- abs(diff(range(at)))/abs(diff(yscale))*vpLocation()$isize["height"]
            acex <- max(0.6, min(vSpaceAvail/vSpaceNeeded, hSpaceAvail/hSpaceNeeded))
        }
        vpTitleAxis <- viewport(x=0.75, width=0.2, yscale=yscale, just=0)
        pushViewport(vpTitleAxis)
        suppressWarnings(grid.yaxis(gp=gpar(col=acol, cex=acex), at=at, label=labs))
        grid.polygon(x=c(0,0,1), y=c(ylim[1],ylim[2],ylim[2]), default.units="native", gp=gpar(col=acol, fill=acol))
        popViewport(2)
    }
    if(.dpOrDefault(GdObject, ".__isCropped", FALSE)){
        vspacing <- as.numeric(convertHeight(unit(2, "points"), "npc"))
        vsize <- as.numeric(convertHeight(unit(4, "points"), "npc"))
        pushViewport(viewport(height=vsize, y=vspacing, just=c(0.5, 0)))
        .moreInd(direction="down", lwd=2)
        popViewport(1)
    }
})


setMethod("drawAxis", signature(GdObject="AlignedReadTrack"), function(GdObject, from, to, subset=TRUE) {
    detail <- match.arg(.dpOrDefault(GdObject, "detail", "coverage"), c("coverage", "reads"))
    if(detail!="coverage") return(NULL) else {
        if(subset)
            GdObject <- subset(GdObject, from=from, to=to)
        cov <- coverage(GdObject, strand="*")
        val <- runValue(coverage(GdObject, strand="*"))
        ## We have to figure out the data range, taking transformation into account
        ylim <- .dpOrDefault(GdObject, "ylim")
        if(is.null(ylim))
        {
            if(!length(val))
                ylim=c(0,1) else{
                    ylim <- c(0, range(val, finite=TRUE, na.rm=TRUE)[2])
                    trans <- .dpOrDefault(GdObject, "transformation")[[1]]
                    if(!is.null(trans))
                        ylim <- c(0, trans(ylim[2]))
                }
        }
        for(s in c("+", "-"))
        {
            pushViewport(viewport(height=0.5, y=ifelse(s=="-", 0, 0.5), just=c("center", "bottom")))
            dummy <- DataTrack(start=rep(mean(c(from, to)),2), end=rep(mean(c(from, to)),2), data=ylim,
                              genome=genome(GdObject), chromosome=chromosome(GdObject))
            oldDp <- displayPars(GdObject, hideInternal=FALSE)
            oldDp[["ylim"]] <- if(s=="+") ylim else rev(ylim)
            displayPars(dummy) <- oldDp
            drawAxis(dummy, from=from, to=to)
            popViewport(1)
        }
    }
})

## Draw a grid in the background of a GdObject. For some subclasses this is meaningless, and the default function will
## return NULL without plotting anything.
setMethod("drawGrid", signature(GdObject="GdObject"), function(GdObject, ...) return(NULL))
setMethod("drawGrid", signature(GdObject="NumericTrack"), function(GdObject, from, to){
    if(.dpOrDefault(GdObject, "grid", FALSE)) {
        vals <- score(GdObject)
        ylim <- .dpOrDefault(GdObject, "ylim", range(vals, na.rm=TRUE, finite=TRUE))
        if(diff(ylim))
        {
            pushViewport(dataViewport(xData=c(from, to), yData=ylim, extension=c(0, 0.1), clip=TRUE))
            panel.grid(h=.dpOrDefault(GdObject, "h", -1), 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)
        }
    }})
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)
    }})
setMethod("drawGrid", signature(GdObject="AlignedReadTrack"), function(GdObject, from, to) {
    detail <- match.arg(.dpOrDefault(GdObject, "detail", "coverage"), c("coverage", "reads"))
    if(detail=="coverage"){
        GdObject <- subset(GdObject, from=from, to=to)
        ## We have to figure out the data range, taking transformation into account
        ylim <- .dpOrDefault(GdObject, "ylim")
        if (is.null(ylim)) {
            maxs <- sapply(c("+", "-"), function(s) {
                cvr <- coverage(GdObject, strand=s)
                                    if (length(cvr)) max(cvr, na.rm=TRUE, finite=TRUE) else 0L
            })
            y.max <- max(maxs, na.rm=TRUE, finite=TRUE)
            ylim <- c(0, if (y.max == 0) 1 else y.max)
            trans <- .dpOrDefault(GdObject, "transformation")[[1]]
            if (!is.null(trans))
                ylim <- c(0, trans(ylim[2]))
        }
        for(s in c("+", "-")) {
            pushViewport(viewport(height=0.5, y=ifelse(s=="-", 0, 0.5), just=c("center", "bottom")))
            dummy <- DataTrack(start=rep(mean(c(from, to)),2), end=rep(mean(c(from, to)),2), data=ylim,
                                                   genome=genome(GdObject), chromosome=chromosome(GdObject))
            oldDp <- displayPars(GdObject, hideInternal=FALSE)
                                oldDp[["ylim"]] <- if(s=="+") ylim else rev(ylim)
            displayPars(dummy) <- oldDp
            drawGrid(dummy, from=from, to=to)
            popViewport(1)
        }
    }
    return(NULL)
})

setMethod("drawGrid", signature(GdObject="AlignmentsTrack"), function(GdObject, from, to){
    if(.dpOrDefault(GdObject, "grid", FALSE)) {
        yvals <- values(GdObject)
        ylim <- .dpOrDefault(GdObject, "ylim", if(!is.null(yvals) && length(yvals))
                             range(yvals, na.rm=TRUE, finite=TRUE) else c(-1,1))
        if(diff(ylim)==0)
            ylim <- ylim+c(-1,1)
        yscale <-c(if(is.null(.dpOrDefault(GdObject, "transformation"))) 0 else min(ylim), max(ylim) + diff(range(ylim)) * 0.05)
        covHeight <- .dpOrDefault(GdObject, ".__coverageHeight", 0)
        pushViewport(viewport(y=1-covHeight["npc"], height=covHeight["npc"], just=c(0.5, 0), yscale=yscale, clip=TRUE))
        panel.grid(v=0, h=.dpOrDefault(GdObject, "h", -1),
                   col=.dpOrDefault(GdObject, "col.grid", "#e6e6e6"), lty=.dpOrDefault(GdObject, "lty.grid", 1),
                   lwd=.dpOrDefault(GdObject, "lwd.grid", 1))
        popViewport(1)
    }
})
##----------------------------------------------------------------------------------------------------------------------------



##----------------------------------------------------------------------------------------------------------------------------
## All the drawGD methods should support two modes, triggered by the boolean argument 'prepare':
##    In prepare mode: nothing is plotted but the object is prepared for plotting bases on the available space. The return
##       value of the method in this mode should always be the updated object. If nothing needs to be prepared, i.e., if the
##       plotting is independent from the available space, simply return the original object
##    In plotting mode: the object is plotted. Return value is the object with optional HTML image map information
##    added to the imageMap slot
## Since subsetting can be potentially expensive when the data are large we want to minimize this operation. Essentially it
## should be done only once before any other plotting or computation starts, hence we expect the GdObject in the drawGD
## methods to already be trimmed to the correct size
##----------------------------------------------------------------------------------------------------------------------------

##----------------------------------------------------------------------------------------------------------------------------
## The default method for all StackedTrack types which always should be called (this has to be done explicitely using
## callNextMethod)
##----------------------------------------------------------------------------------------------------------------------------
## Although the stacking type is not stored as a displayParameter we still want to check whether it is
## included there and set the actual stacking of the object accordingly
setMethod("drawGD", signature("StackedTrack"), function(GdObject, ...){
    debug <- .dpOrDefault(GdObject, "debug", FALSE)
    if((is.logical(debug) && debug))
        browser()
    st <- .dpOrDefault(GdObject, "stacking")
    if(!is.null(st))
        stacking(GdObject) <- st
    return(invisible(GdObject))
})
##----------------------------------------------------------------------------------------------------------------------------


##----------------------------------------------------------------------------------------------------------------------------
## The plotting method for a CustomTrack actually just calls the user-defined plotting function
##----------------------------------------------------------------------------------------------------------------------------
## Although the stacking type is not stored as a displayParameter we still want to check whether it is
## included there and set the actual stacking of the object accordingly
setMethod("drawGD", signature("CustomTrack"), function(GdObject, minBase, maxBase, prepare=FALSE, ...){
    rev <- .dpOrDefault(GdObject, "reverseStrand", FALSE)
    xscale <- if(!rev) c(minBase, maxBase) else c(maxBase, minBase)
    pushViewport(viewport(xscale=xscale, clip=TRUE))
    tmp <- GdObject@plottingFunction(GdObject, prepare=prepare)
    if(!is(tmp, "CustomTrack")){
        warning("The plotting function of a CustomTrack has to return the input object. Using the original CustomTrack object now.")
    }else{
        GdObject <- tmp
    }
    popViewport(1)
    return(invisible(GdObject))
})
##----------------------------------------------------------------------------------------------------------------------------


##----------------------------------------------------------------------------------------------------------------------------
## The method for OverlayTracks simply delegates to the methods for each of the elements in trackList
##----------------------------------------------------------------------------------------------------------------------------
setMethod("drawGD", signature("OverlayTrack"), function(GdObject, ...){
    GdObject@trackList <- lapply(GdObject@trackList, drawGD, ...)
    return(invisible(GdObject))
})
##----------------------------------------------------------------------------------------------------------------------------


##----------------------------------------------------------------------------------------------------------------------------
## 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
.barsAndLabels <- function(GdObject){
    bins <- stacks(GdObject)
    stacks <- max(bins)
    res <- .pxResolution(coord="x")
    gp <- group(GdObject)
    grpSplit <- split(range(GdObject), gp)
    grpRanges <- unlist(range(grpSplit))
    needBar <- sapply(grpSplit, length)>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 substract 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 <- sapply(split(strand(GdObject), gp), function(x){
        tmp <- unique(x)
        if(length(tmp)>1) "*" else tmp
    })
    yloc <- sapply(split((stacks-bins)+1, gp), function(x) unique(x))+0.5
    color <- if(length(grep("__desatCol", values(GdObject)$feature[1])))
        .dpOrDefault(GdObject, "fill", .DEFAULT_FILL_COL) else sapply(split(.getBiotypeColor(GdObject), gp), head, 1)
    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(sapply(split(.getAnn(GdObject, "gdensity"), gp), head, 1))
            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)))
                   cy <- yloc[indx] + 0.5
                   labs <- labs[indx]
                   algn <- c("center", "top")
               },
               "below"={
                   cx <- start(grpRangesCut) + width(grpRangesCut)/2
                   indx <- which(seq_len(length(grpRanges)) %in% queryHits(findOverlaps(grpRanges, grpRangesCut)))
                   cy <- yloc[indx] - 0.5
                   labs <- labs[indx]
                   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))
}

## The actual drawing method
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) && !is.na(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[,1:4, drop=FALSE],
                      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 elipses
        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(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)
    ## Finaly 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 <- sapply(restCols, function(x){
            tmp <- as.character(box[,x])
            names(tmp) <- rownames(coords)
            tmp}, simplify=FALSE)
        tags$title <- identifier(GdObject)
        ImageMap(coords=coords, tags=tags) } else NULL
    imageMap(GdObject) <- im
    return(invisible(GdObject))
})


## For a GeneRegionTrack we just set the showExonId alias and then call the AnnotationTrack method
setMethod("drawGD", signature("GeneRegionTrack"), function(GdObject,  ...){
    displayPars(GdObject) <- list(showFeatureId=as.vector(.dpOrDefault(GdObject, "showExonId")))
    GdObject <- callNextMethod()
    return(invisible(GdObject))
})

## The actual drawing method
setMethod("drawGD", signature("AlignmentsTrack"), function(GdObject, minBase, maxBase, prepare=FALSE, subset=TRUE, ...){
    debug <- .dpOrDefault(GdObject, "debug", FALSE)
    imageMap(GdObject) <- NULL
    if(!length(GdObject))
        return(invisible(GdObject))
    rev <- .dpOrDefault(GdObject, "reverseStrand", FALSE)
    revs <- !.dpOrDefault(GdObject, "reverseStacking", FALSE)
    vSpacing <- as.numeric(convertHeight(unit(3, "points"), "npc"))
    type <- match.arg(.dpOrDefault(GdObject, "type", .ALIGNMENT_TYPES), .ALIGNMENT_TYPES, several.ok=TRUE)
    ylim <- .dpOrDefault(GdObject, "ylim")
    ## 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 compute the coverage vector which will be needed for the axis
    if(prepare){
        if((is.logical(debug) && debug) || "prepare" %in% debug)
            browser()
        GdObject <- callNextMethod(GdObject, ...)
        ## The mismatched bases need to be extracted from the read sequences and the reference sequence
        if(.dpOrDefault(GdObject, "showMismatches", TRUE) && !is.null(GdObject@referenceSequence)){
            mm <- .findMismatches(GdObject)
            if(nrow(mm))
                displayPars(GdObject) <- list(".__mismatches"=mm)
        }
        ## The coverage calculation and the height of the coverage section
        if("coverage" %in% type){
            covSpace <- as.numeric(convertHeight(unit(5, "points"), "npc"))
            if("pileup" %in% type){
                covHeight <- .dpOrDefault(GdObject, "coverageHeight", 0.1)
                if(covHeight > 0 && covHeight < 1)
                    covHeight <- as.numeric(convertHeight(unit(covHeight, "npc"), "points"))
                covHeight <- max(.dpOrDefault(GdObject, "minCoverageHeight", 50), covHeight)
                covHeight <- c(points=covHeight, npc=as.numeric(convertHeight(unit(covHeight, "points"), "npc")))
            }else if ("sashimi" %in% type){
                covHeight <- c(npc=0.5-(covSpace*2), points=1)
            }else{
                covHeight <- c(npc=1-(covSpace*2), points=1)
            }
            coverage <- coverage(range(GdObject), width=maxBase)
            ## apply data transformation if one is set up
            trans <- .dpOrDefault(GdObject, "transformation")
            if(is.list(trans))
                trans <- trans[[1]]
            if(!is.null(trans)){
                if(!is.function(trans) || length(formals(trans))!=1L)
                    stop("Display parameter 'transformation' must be a function with a single argument")
                test <- trans(coverage)
                if(!is(test, "Rle") || length(test) != length(coverage))
                    stop("The function in display parameter 'transformation' results in invalid output.\n",
                         "It has to return a numeric matrix with the same dimensions as the input data.")
                coverage <- test
            }
            displayPars(GdObject) <- list(".__coverage"=coverage,
                                          ".__coverageHeight"=covHeight,
                                          ".__coverageSpace"=covSpace)
        } else {
            covHeight <- c(npc=0, points=0)
            covSpace <- 0
        }
        ## The sashimi calculation and the height of the sashimi section
        if("sashimi" %in% type){
            sashSpace <- as.numeric(convertHeight(unit(5, "points"), "npc"))
            if("pileup" %in% type) {
                sashHeight <- .dpOrDefault(GdObject, "sashimiHeight", 0.1)
                if(sashHeight > 0 && sashHeight < 1)
                    sashHeight <- as.numeric(convertHeight(unit(sashHeight, "npc"), "points"))
                sashHeight <- max(.dpOrDefault(GdObject, "minSashimiHeight", 50), sashHeight)
                sashHeight <- c(points=sashHeight, npc=as.numeric(convertHeight(unit(sashHeight, "points"), "npc")))
            }else if ("coverage" %in% type){
                sashHeight <- c(npc=0.5-(sashSpace*2), points=1)
            }else{
                sashHeight <- c(npc=1-(sashSpace*2), points=1)
            }
            sashScore <- .dpOrDefault(GdObject, "sashimiScore", 1L)
            sashLwdMax <- .dpOrDefault(GdObject, "lwd.sashimiMax", 10)
            sashStrand <- .dpOrDefault(GdObject, "sashimiStrand", "*")
            sashFilter <- .dpOrDefault(GdObject, "sashimiFilter", NULL)
            sashFilterTolerance <- .dpOrDefault(GdObject, "sashimiFilterTolerance", 0L)
            sashNumbers <- .dpOrDefault(GdObject, "sashimiNumbers", FALSE)
            sash <- .dpOrDefault(GdObject, "sashimiJunctions", NULL)
            if (is.null(sash)) {
                sash <- .create.summarizedJunctions.for.sashimi.junctions(ranges(GdObject))
            } else {
                if (class(sash)!="GRanges")
                    stop("\"sashimiJunctions\" object must be of \"GRanges\" class!")
                sashMcolName <- if (sashStrand=="+") "plus_score" else if (sashStrand=="-") "minus_score" else "score"
                if (sum(colnames(mcols(sash))==sashMcolName) != 1)
                    stop(sprintf("\"mcols\" of \"sashimiJunctions\" object must contain column named \"%s\",\n which matches the specified (%s) \"sashimiStrand\"!", sashMcolName, sashStrand))
            }
            sash <- .convert.summarizedJunctions.to.sashimi.junctions(juns=sash,
                                                                      score=sashScore,
                                                                      lwd.max=sashLwdMax,
                                                                      strand=sashStrand,
                                                                      filter=sashFilter,
                                                                      filterTolerance=sashFilterTolerance)
            displayPars(GdObject) <- list(".__sashimi"=sash,
                                          ".__sashimiHeight"=sashHeight,
                                          ".__sashimiSpace"=sashSpace,
                                          ".__sashimiNumbers"=sashNumbers)
        } else {
            sashHeight <- c(npc=0, points=0)
            sashSpace <- 0
        }
        if("pileup" %in%