#' @include DataTrack-class.R
#' @include ReferenceTrack-class.R
NULL
## AnnotationTrack Class -----------------------------------------------------
#' AnnotationTrack class and methods
#'
#' A fairly generic track object for arbitrary genomic range annotations, with
#' the option of grouped track items. The extended `DetailsAnnotationTrack`
#' provides a more flexible interface to add user-defined custom
#' information for each range.
#'
#' @name AnnotationTrack-class
#'
#' @template GdObject-class_slot
#' @template RangeTrack-class_slot
#' @template StackedTrack_slot
#' @template AnnotationTrack_slot
#'
#' @template AnnotationTrack_param
#'
#' @author Florian Hahne, Arne Mueller
#'
#' @inherit GdObject-class seealso
#'
#' @examples
#' ## An empty object
#' AnnotationTrack()
#'
#' ## Construct from individual arguments
#' st <- c(2000000, 2070000, 2100000, 2160000)
#' ed <- c(2050000, 2130000, 2150000, 2170000)
#' str <- c("-", "+", "-", "-")
#' gr <- c("Group1", "Group2", "Group1", "Group3")
#'
#' annTrack <- AnnotationTrack(
#' start = st, end = ed, strand = str, chromosome = 7,
#' genome = "hg19", feature = "test", group = gr,
#' id = paste("annTrack item", 1:4),
#' name = "generic annotation", stacking = "squish"
#' )
#'
#' ## Or from a data.frame
#' df <- data.frame(
#' start = st, end = ed, strand = str, id = paste("annTrack item", 1:4),
#' feature = "test", group = gr
#' )
#' annTrack <- AnnotationTrack(
#' range = df, genome = "hg19", chromosome = 7,
#' name = "generic annotation", stacking = "squish"
#' )
#'
#' ## Or from a GRanges object
#' gr <- GRanges(
#' seqnames = "chr7", range = IRanges(start = df$start, end = df$end),
#' strand = str
#' )
#' genome(gr) <- "hg19"
#' mcols(gr) <- df[, -(1:3)]
#' annTrack <- AnnotationTrack(
#' range = gr, name = "generic annotation",
#' stacking = "squish"
#' )
#'
#' ## Finally from a GRangesList
#' grl <- split(gr, values(gr)$group)
#' AnnotationTrack(grl)
#' \dontshow{
#' ## For some annoying reason the postscript device does not know about
#' ## the sans font
#' if (!interactive()) {
#' font <- ps.options()$family
#' displayPars(annTrack) <- list(fontfamily = font, fontfamily.title = font)
#' }
#' }
#'
#' ## Plotting
#' plotTracks(annTrack)
#'
#' ## Track names
#' names(annTrack)
#' names(annTrack) <- "foo"
#' plotTracks(annTrack)
#'
#' ## Subsetting and splitting
#' subTrack <- subset(annTrack, to = 2155000)
#' length(subTrack)
#' subTrack[1:2]
#' split(annTrack, c(1, 2, 1, 2))
#'
#' ## Accessors
#' start(annTrack)
#' end(annTrack)
#' width(annTrack)
#' position(annTrack)
#' width(subTrack) <- width(subTrack) + 1000
#'
#' strand(annTrack)
#' strand(subTrack) <- "-"
#'
#' chromosome(annTrack)
#' chromosome(subTrack) <- "chrX"
#'
#' genome(annTrack)
#' genome(subTrack) <- "mm9"
#'
#' range(annTrack)
#' ranges(annTrack)
#'
#' ## Annotation
#' identifier(annTrack)
#' identifier(annTrack, "lowest")
#' identifier(subTrack) <- "bar"
#'
#' feature(annTrack)
#' feature(subTrack) <- "foo"
#'
#' values(annTrack)
#'
#' ## Grouping
#' group(annTrack)
#' group(subTrack) <- "Group 1"
#' chromosome(subTrack) <- "chr7"
#' plotTracks(subTrack)
#'
#' ## Stacking
#' stacking(annTrack)
#' stacking(annTrack) <- "dense"
#' plotTracks(annTrack)
#'
#' ## coercion
#' as(annTrack, "data.frame")
#' as(annTrack, "UCSCData")
#'
#' ## HTML image map
#' coords(annTrack)
#' tags(annTrack)
#' annTrack <- plotTracks(annTrack)$foo
#' coords(annTrack)
#' tags(annTrack)
#'
#' ## DetailsAnnotationTrack
#' library(lattice) # need to use grid grapics
#'
#' ## generate two random distributions per row (probe/feature)
#' ## the difference between the distributions increases from probe 1 to 4
#' m <- matrix(c(rgamma(400, 1)), ncol = 100)
#' m[, 51:100] <- m[, 51:100] + 0:3
#' ## rownames must be accessible by AnnotationTrack element identifier
#' rownames(m) <- identifier(annTrack, "lowest")
#'
#' ## create a lattice density plot for the values (signals) of the two groups
#' ## as the chart must be placed into a pre-set grid view port we have to use
#' ## print without calling plot.new! Note, use a common prefix for all lattice.
#' ## Avoid wasting space by removing y-axis decorations.
#'
#' ## Note, in this example 'm' will be found in the environment the 'details'
#' ## function is defined in. To avoid overwriting 'm' you should use a closure
#' ## or environment to access 'm'.
#' details <- function(identifier, ...) {
#' d <- data.frame(signal = m[identifier, ], group = rep(c("grp1", "grp2"), each = 50))
#' print(densityplot(~signal,
#' group = group, data = d, main = identifier,
#' scales = list(draw = FALSE, x = list(draw = TRUE)), ylab = "", xlab = "",
#' ), newpage = FALSE, prefix = "plot")
#' }
#'
#' deTrack <- AnnotationTrack(
#' range = gr, genome = "hg19", chromosome = 7,
#' name = "generic annotation with details per entry", stacking = "squish",
#' fun = details, details.ratio = 1
#' )
#'
#' plotTracks(deTrack)
#'
#' set.seed(1234)
#' deTrack <- AnnotationTrack(
#' range = gr, genome = "hg19", chromosome = 7,
#' name = "generic annotation with details per entry",
#' stacking = "squish", fun = details,
#' details.ratio = 1, selectFun = function(...) {
#' sample(c(FALSE, TRUE), 1)
#' }
#' )
#'
#' plotTracks(deTrack)
#' @importClassesFrom rtracklayer UCSCData
#'
#' @exportClass AnnotationTrack
setClass("AnnotationTrack",
contains = "StackedTrack",
prototype = prototype(
columns = c("feature", "group", "id"),
stacking = "squish",
name = "AnnotationTrack",
dp = DisplayPars(
arrowHeadWidth = 30,
arrowHeadMaxWidth = 40,
cex.group = 0.6,
cex = 1,
col.line = "darkgray",
col = "transparent",
featureAnnotation = NULL,
fill = "lightblue",
fontfamily.group = "sans",
fontcolor.group = .DEFAULT_SHADED_COL,
fontcolor.item = "white",
fontface.group = 2,
fontsize.group = 12,
groupAnnotation = NULL,
just.group = "left",
lex = 1,
lineheight = 1,
lty = "solid",
lwd = 1,
mergeGroups = FALSE,
min.height = 3,
min.width = 1,
rotation = 0,
rotation.group = 0,
rotation.item = 0,
shape = "arrow",
showFeatureId = FALSE,
showId = FALSE,
showOverplotting = FALSE,
size = 1
)
)
)
## Initialize ----------------------------------------------------------------
## Essentially we just check for the correct GRanges columns here
#' @describeIn AnnotationTrack-class Show method.
#' @export
setMethod("initialize", "AnnotationTrack", function(.Object, ...) {
if (is.null(list(...)$range) && is.null(list(...)$genome) && is.null(list(...)$chromosome)) {
return(.Object)
}
## the display parameter defaults
.makeParMapping()
.Object <- .updatePars(.Object, "AnnotationTrack")
range <- list(...)$range
if (!is.null(range) && length(.Object)) {
if (!all(.Object@columns %in% colnames(values(range)))) {
stop(
"Problem initializing AnnotationTrack, need the following columns:",
toString(.Object@columns, collapse = ", ")
)
}
grp <- if (is(.Object, "GeneRegionTrack")) values(range)$transcript else values(range)$group
if (any(vapply(split(as.character(strand(range)), grp), function(x) length(unique(x)), numeric(1)) != 1)) {
stop("Grouped elments of a RangeTrack can not be on opposing strands")
}
}
.Object <- callNextMethod()
return(.Object)
})
## The file-based version of the AnnotationTrack class. This will mainly provide a means to dispatch to
## a special 'subset' method which should stream the necessary data from disk.
#' @describeIn AnnotationTrack-class The file-based version of the `AnnotationTrack-class`.
#' @exportClass ReferenceAnnotationTrack
setClass("ReferenceAnnotationTrack", contains = c("AnnotationTrack", "ReferenceTrack"))
## Initialize ----------------------------------------------------------------
## This just needs to set the appropriate slots that are being inherited from ReferenceTrack because the
## multiple inheritance has some strange features with regards to method selection
#' @describeIn AnnotationTrack-class Initialize.
#' @export
setMethod("initialize", "ReferenceAnnotationTrack", function(.Object, stream, reference, mapping = list(),
args = list(), defaults = list(), ...) {
.Object <- selectMethod("initialize", "ReferenceTrack")(.Object = .Object, reference = reference, stream = stream,
mapping = mapping, args = args, defaults = defaults)
.Object <- callNextMethod(.Object, ...)
return(.Object)
})
## Constructor ---------------------------------------------------------------
## Constructor. The following arguments are supported:
## o range: a data.frame or a GRanges object containing the information
## about the track items. If a data.frame, it needs to be coerceable
## to a GRanges object, i.e., it needs at least the mandatory 'start', 'stop' and
## 'strand' columns. Additional optional columns are:
## - feature: the type of the item. Can be mapped to colors via the DisplayPars.
## - group: a grouping factor to connect track items.
## - id: a unique identifier for a feature. This will be plotted if showId==TRUE
## Note that internally we use the value of ID as the seqnames slot in the
## internal GRanges object. Defaults for all missing columns are generated.
## Instead of using the 'range' parameter, all these values can also be passed as
## individual vectors, in which case they need to be of similar length.
## o start, end, width: numeric vectors of the item start and end coordinates
## o strand: the strand information may be provided in the form '+' for
## the Watson strand, '-' for the Crick strand or '*' for any of the two.
## o feature, group, id: individual vectors of equal length as described above.
## o genome, chromosome: the reference genome and active chromosome for the track.
## o stacking: character controlling the stacking of overlapping items. One in 'hide',
## 'dense', 'squish', 'pack' or 'full'.
## o name: the name of the track. This will be used for the title panel.
## All additional items in ... are being treated as further DisplayParameters
#' @describeIn AnnotationTrack-class Constructor function for
#' `AnnotationTrack-class`
#' @return
#' The return value of the constructor function is a new object of class
#' \code{AnnotationTrack} or of class \code{DetailsAnnotationTrack}, depending
#' on the constructor arguments. Typically the user will not have to be
#' troubled with this distinction and can rely on the constructor to make the
#' right choice.
#'
#' @section Objects from the class:
#'
#' Objects can be created using the constructor function
#' \code{AnnotationTrack}.
#' @export
AnnotationTrack <- function(range = NULL, start = NULL, end = NULL, width = NULL, feature, group, id, strand, chromosome,
genome, stacking = "squish", name = "AnnotationTrack", fun, selectFun, importFunction,
stream = FALSE, ...) {
## Some defaults
covars <- .getCovars(range)
isStream <- FALSE
if (!is.character(range)) {
n <- max(c(length(start), length(end), length(width)), nrow(covars))
if (is.null(covars[["feature"]]) && missing(feature)) {
feature <- rep("unknown", n)
}
if (is.null(covars[["id"]]) && missing(id)) {
id <- make.unique(rep(if (!is.null(feature)) as.character(feature) else covars[["feature"]], n)[seq_len(n)])
}
if (is.null(covars[["group"]]) && missing(group)) {
group <- seq_len(n)
}
}
## Build a GRanges object from the inputs
.missingToNull(c("feature", "group", "id", "strand", "chromosome", "importFunction", "genome"))
args <- list(feature = feature, group = group, id = id, strand = strand, chromosome = chromosome, genome = genome)
defs <- list(feature = "unknown", group = "unknown", id = "unknown", strand = "*", density = 1, chromosome = "chrNA", genome = NA)
range <- .buildRange(
range = range, groupId = "group", start = start, end = end, width = width,
args = args, defaults = defs, chromosome = chromosome, trackType = "AnnotationTrack",
importFun = importFunction, stream = stream
)
if (is.list(range)) {
isStream <- TRUE
slist <- range
range <- GRanges()
}
## Pipes have a special meaning for merged groups, so we can't have them in the initial group vector
mcols(range)[["group"]] <- gsub("|", "", mcols(range)[["group"]], fixed = TRUE)
## If no chromosome was explicitely asked for we just take the first one in the GRanges object
if (missing(chromosome) || is.null(chromosome)) {
chromosome <- if (length(range) > 0) .chrName(as.character(seqnames(range)[1])) else "chrNA"
}
## And finally the object instantiation, we have to distinguish between DetailsAnnotationTracks and normal ones
genome <- .getGenomeFromGRange(range, ifelse(is.null(genome), character(), genome[1]))
if (missing(fun)) {
if (!isStream) {
return(new("AnnotationTrack",
chromosome = as.character(chromosome[1]), range = range,
name = name, genome = genome, stacking = stacking, ...
))
} else {
## A bit hackish but for some functions we may want to know which track type we need but at the
## same time we do not want to enforce this as an additional argument
e <- new.env()
e[["._trackType"]] <- "AnnotationTrack"
environment(slist[["stream"]]) <- e
return(new("ReferenceAnnotationTrack",
chromosome = as.character(chromosome[1]), range = range,
name = name, genome = genome, stacking = stacking, stream = slist[["stream"]], reference = slist[["reference"]],
mapping = slist[["mapping"]], args = args, defaults = defs, ...
))
}
} else {
if (!is.function(fun)) {
stop("'fun' must be a function")
}
if (missing(selectFun)) {
selectFun <- function(...) {
return(TRUE)
}
}
if (!is.function(selectFun)) {
stop("'selectFun' must be a function")
}
return(new("DetailsAnnotationTrack",
chromosome = as.character(chromosome[1]), range = range,
name = name, genome = genome, stacking = stacking, fun = fun, selectFun = selectFun, ...
))
}
}
#' @describeIn AnnotationTrack-class directly extends `AnnotationTrack.`
#' @exportClass DetailsAnnotationTrack
setClass("DetailsAnnotationTrack",
contains = "AnnotationTrack",
representation = representation(fun = "function", selectFun = "function"),
prototype = prototype(
fun = function(...) {},
selectFun = function(...) {
return(TRUE)
},
dp = DisplayPars(
details.minWidth = 100,
details.ratio = Inf,
details.size = 0.5,
detailsBorder.col = "darkgray",
detailsBorder.fill = "transparent",
detailsBorder.lty = "solid",
detailsBorder.lwd = 1,
detailsConnector.cex = 1,
detailsConnector.col = "darkgray",
detailsConnector.lty = "dashed",
detailsConnector.lwd = 1,
detailsConnector.pch = 20,
detailsFunArgs = list(),
groupDetails = FALSE
)
)
)
#' @describeIn AnnotationTrack-class Constructor function for
#' `DetailsAnnotationTrack-class`
#'
#' The `DetailsAnnotationTrack` class directly extends `AnnotationTrack.`
#' The purpose of this track type is to add an arbitrarily detailed plot
#' section (typically consisting of additional quantitative data) for each
#' range element of an `AnnotationTrack.` This allows a locus wide view of
#' annotation elements together with any kind of details per feature or element
#' that may for instance provide insight on how some complex quantitative
#' measurements change according to their position in a locus. If the
#' quantitative data is too complex for a `DataTrack` e.g. because it requires
#' extra space or a trellis-like representation, a `DetailsAnnotationTrack` can
#' be used instead. Example: An `AnnotationTrack` shows the positions of a
#' number of probes from a microarray, and you want a histogram of the signal
#' intensity distribution derived from all samples at each of these probe
#' location. Another example usage would be to show for each element of an
#' `AnnotationTrack` an xy-plot of the signal against some clinical measurement
#' such as blood pressure. The limitation for applications of this type of
#' track is basically only the available space of the device you are
#' plotting to.
#'
#' This flexibility is possible by utilizing a simple function model
#' to perform all the detailed plotting. The functionality of this plotting
#' function fun is totally up to the user, and the function environment is
#' prepared in a way that all necessary information about the plotted
#' annotation feature is available. To restrict the details section to only
#' selected number of annotation features one can supply another function
#' `selectFun`, which decides for each feature separately whether details are
#' available or not. Finally, an arbitrary number of additional arguments can
#' be passed on to these two function by means of the `detailsFunArgs` display
#' parameter. This is expected to be a named list, and all list elements are
#' passed along to the plotting function fun and to the selector function
#' `selectFun` as additional named arguments. Please note that some argument
#' names like `start`, `end` or `identifier` are reserved and can not be used
#' in the `detailsFunArgs` list. For examples of plotting functions,
#' see the 'Examples' section.
#' @export
DetailsAnnotationTrack <- function(...) AnnotationTrack(...)
## Initialize ----------------------------------------------------------------
#' @describeIn AnnotationTrack-class Initialize.
#' @export
setMethod("initialize", "DetailsAnnotationTrack", function(.Object, fun, selectFun, ...) {
## the diplay parameter defaults
.Object <- .updatePars(.Object, "DetailsAnnotationTrack")
.makeParMapping()
.Object@fun <- fun
.Object@selectFun <- selectFun
.Object <- callNextMethod(.Object, ...)
return(.Object)
})
## General accessors ---------------------------------------------------------
## Annotation Accessors ------------------------------------------------------
#' @describeIn AnnotationTrack-class extract the group membership for all track items.
#' @export
setMethod("group", "AnnotationTrack", function(GdObject) .getAnn(GdObject, "group"))
#' @describeIn AnnotationTrack-class replace the grouping information for track items.
#' The replacement value must be a factor of appropriate length or another
#' vector that can be coerced into such.
#' @export
setReplaceMethod("group", signature("AnnotationTrack", "character"), function(GdObject, value) .setAnn(GdObject, value, "group"))
## Context-dependent meta-accessors to the identifier data of an AnnotationTrack and a GeneRegionTrack. For the former, those will
## be the content of the 'id', the 'group' or the 'feature' metadata column, for the latter, one in
## the selection of 'symbol', 'gene', 'transcript', 'feature' or 'exon'. For historical reasons, logical values are also allowed, where TRUE
## is equivalent to 'symbol' and FALSE is equivalent to 'gene'. If not provided explicitely in the 'type' argument, the identifier
## type is inferred from the 'transcriptAnnotation' display parameter (a.k.a. 'geneSymbols') for a GeneRegionTrack, and from the
## 'groupAnnotation' parameter for an AnnotationTrack. The special value 'lowest' should work across all classes and provided the most
## detailed annotation ('id' for AnnotationTrack and 'exon' for GeneRegionTrack).
#' @describeIn AnnotationTrack-class return track item identifiers.
#' Depending on the setting of the optional argument lowest, these are either
#' the group identifiers or the individual item identifiers.
#' @export
setMethod("identifier", "AnnotationTrack", function(GdObject, type = .dpOrDefault(GdObject, "groupAnnotation", "group")) {
if (is.null(type)) {
type <- "group"
}
id <- switch(as.character(type),
"group" = group(GdObject),
"id" = .getAnn(GdObject, "id"),
"feature" = feature(GdObject),
"lowest" = .getAnn(GdObject, "id"),
group(GdObject)
)
id[is.na(id)] <- "NA"
return(id)
})
#' @describeIn AnnotationTrack-class Set the track item identifiers.
#' The replacement value has to be a character vector of appropriate length.
#' This always replaces the group-level identifiers, so essentially it is
#' similar to `groups<-`.
#' @export
setReplaceMethod("identifier", c("AnnotationTrack", "character"), function(GdObject, value) {
type <- .dpOrDefault(GdObject, "groupAnnotation", "group")
switch(as.character(type),
"group" = group(GdObject) <- value,
"id" = ranges(GdObject)$id <- value,
"feature" = feature(GdObject) <- value,
"lowest" = ranges(GdObject)$id <- value,
group(GdObject) <- value
)
return(GdObject)
})
## Stacking ------------------------------------------------------------------
#' @describeIn AnnotationTrack-class Recompute the stacks based on the available
#' space and on the object's track items and stacking settings.
#' @export
setMethod("setStacks", "AnnotationTrack", function(GdObject, recomputeRanges = TRUE) {
if (!.needsStacking(GdObject) || length(GdObject) == 0) {
## No stacks needed, so there is just a single bin
bins <- rep(1, length(GdObject))
} else {
gp <- group(GdObject)
if (!is.factor(gp)) {
gp <- factor(gp, levels=unique(gp))
}
needsGrp <- any(duplicated(gp))
lranges <- .dpOrDefault(GdObject, ".__groupRanges")
gpt <- if (needsGrp) table(gp) else rep(1, length(GdObject))
if (recomputeRanges || is.null(lranges) || length(lranges) != length(gpt)) {
GdObject <- .computeGroupRange(GdObject)
lranges <- .dpOrDefault(GdObject, ".__groupRanges")
}
uid <- if (.transcriptsAreCollapsed(GdObject)) {
sprintf("uid%i", seq_along(identifier(GdObject)))
} else {
make.unique(identifier(GdObject, type = "lowest"))
}
bins <- rep(disjointBins(lranges), gpt)
names(bins) <- if (needsGrp) unlist(split(uid, gp)) else uid
bins <- bins[uid]
}
bins <- if (length(bins)) bins else 0
GdObject@stacks <- bins
return(GdObject)
})
## Consolidate ---------------------------------------------------------------
#' @describeIn AnnotationTrack-class Consolidate.
#' Determine whether there is group label annotation or not, and add this
#' information as the internal display parameter `.__hasAnno`. Precompute
#' the grouped ranges together with optional labels in order to determine
#' the correct plotting range later.
#' @keywords internal
setMethod("consolidateTrack", signature(GdObject = "AnnotationTrack"), function(GdObject, hasAxis = FALSE,
hasTitle = .dpOrDefault(GdObject, "showTitle", TRUE),
title.width = NULL, ...) {
GdObject <- callNextMethod(GdObject, ...)
if (length(GdObject)) {
## ids are shown if either set by the showId parameter, or through the use of the transcriptAnnotation or
## groupAnnotation parameters
ids <- identifier(GdObject)
sid <- .dpOrDefault(GdObject, "showId")
ta <- if (is(GdObject, "GeneRegionTrack")) .dpOrDefault(GdObject, "transcriptAnnotation") else .dpOrDefault(GdObject, "groupAnnotation")
sid <- if (is.null(sid)) !is.null(ta) && ta != "none" else sid
hasAnno <- sid && !all(ids == "")
## element ids are shown if either the showFeatureId or showExonId parameters are TRUE, or if featureAnnotation
## or exonAnnotation is set
sfid <- if (is(GdObject, "GeneRegionTrack")) .dpOrDefault(GdObject, "showExonId") else .dpOrDefault(GdObject, "showFeatureId")
fa <- if (is(GdObject, "GeneRegionTrack")) .dpOrDefault(GdObject, "exonAnnotation") else .dpOrDefault(GdObject, "featureAnnotation")
sfid <- if (is.null(sfid)) !is.null(fa) && fa != "none" else sfid
displayPars(GdObject) <- list(".__hasAnno" = hasAnno, showId = sid, showFeatureId = sfid, showExonId = sfid)
GdObject <- .computeGroupRange(GdObject, hasAxis = hasAxis, hasTitle = hasTitle, title.width = title.width)
} else {
displayPars(GdObject) <- list(".__hasAnno" = FALSE, showId = FALSE, showFeatureId = FALSE, showExonId = FALSE)
}
return(GdObject)
})
## Collapse -----------------------------------------------------------------
## A slightly quicker function to compute overlaps between two GRanges objects
.myFindOverlaps <- function(gr1, gr2) {
gr1 <- sort(gr1)
gr2 <- sort(gr2)
gr1 <- split(gr1, seqnames(gr1))
gr2 <- split(gr2, seqnames(gr2))
queryHits(findOverlaps(ranges(gr1), ranges(gr2)))
}
## Find all elements in a GRanges object 'grange' with distance smaller than 'minXDist' and merge them along with their additional
## metadata columns. 'elements' is a frequency table of items per group, and it is needed to figure out whether all items of a given
## group have been merged. 'GdObject' is the input track object from which certain information has to be extracted. The output of this
## function is a list with elements
## o range: the updated GRanges object
## o needsRestacking: logical flag indicating whether stacks have to be recomputed
## o split: the original merged annotation (need to work on this, currently not used)
## o merged: logical vector indicating which of the elements in 'range' constitute fully merged groups
.collapseAnnotation <- function(grange, minXDist, elements, GdObject, offset = 0) {
needsRestacking <- TRUE
annoSplit <- merged <- NULL
anno <- as.data.frame(grange)
for (i in colnames(anno)) {
if (is.factor(anno[, i])) {
anno[, i] <- as.character(anno[, i])
}
}
cols <- c("strand", "density", "gdensity", "feature", "id", "start", "end", if (is(GdObject, "GeneRegionTrack")) {
c("gene", "exon", "transcript", "symbol", "rank")
} else {
"group"
})
missing <- which(!cols %in% colnames(anno))
for (i in missing) {
anno[, cols[missing]] <- if (cols[i] == "density") 1 else NA
}
rRed <- if (length(grange) > 1) reduce(grange, min.gapwidth = minXDist, with.revmap = TRUE) else grange
if (length(rRed) < length(grange)) {
## Some of the items have to be merged and we need to make sure that the additional annotation data that comes with it
## is processed in a sane way.
needsRestacking <- TRUE
mapping <- rep(seq_along(rRed$revmap), elementNROWS(rRed$revmap))
## We start by finding the items that have not been reduced
identical <- mapping %in% which(table(mapping) == 1)
newVals <- anno[identical, cols]
## Here we hijack the seqnames column to indicate whether the whole group has been merged
if (nrow(newVals)) {
newVals$seqnames <- elements[as.character(anno[identical, "seqnames"])] == 1
newVals$gdensity <- ifelse(elements[as.character(anno[identical, "seqnames"])] == 1, 1, NA)
}
## Now find out which original items have been merged
grange <- grange[!identical]
rRed <- rRed[-(mapping[identical])]
index <- mapping[!identical]
annoSplit <- split(anno[!identical, ], index)
cid <- function(j) sprintf("[Cluster_%i] ", j + offset)
## FIXME: We could speed this up by running it in C
tmpVals <- lapply(seq_along(annoSplit), function(i) {
x <- annoSplit[[i]]
if (is(GdObject, "GeneRegionTrack")) {
c(
strand = ifelse(length(unique(x[, "strand"])) == 1, as.character(x[1, "strand"]), "*"),
density = sum(as.integer(x[, "density"])),
gdensity = ifelse(is.na(head(x[, "gdensity"], 1)), 1, sum(as.integer(x[, "gdensity"]))),
feature = ifelse(length(unique(x[, "feature"])) == 1, as.character(x[1, "feature"]), "composite"),
id = ifelse(length(unique(x[, "id"])) == 1, as.character(x[1, "id"]), cid(i)),
start = min(x[, "start"]),
end = max(x[, "end"]),
gene = ifelse(length(unique(x[, "gene"])) == 1, as.character(x[1, "gene"]), cid(i)),
exon = ifelse(length(unique(x[, "exon"])) == 1, as.character(x[1, "exon"]), cid(i)),
transcript = ifelse(length(unique(x[, "transcript"])) == 1, as.character(x[1, "transcript"]), cid(i)),
symbol = ifelse(length(unique(x[, "symbol"])) == 1, as.character(x[1, "symbol"]), cid(i)),
rank = min(as.integer(x[, "rank"])), seqnames = as.vector(nrow(x) == elements[x[1, "seqnames"]])
)
} else {
c(
strand = ifelse(length(unique(x[, "strand"])) == 1, as.character(x[1, "strand"]), "*"),
density = sum(as.integer(x[, "density"])),
gdensity = ifelse(is.na(head(x[, "gdensity"], 1)), 1, sum(as.integer(x[, "gdensity"]))),
feature = ifelse(length(unique(x[, "feature"])) == 1, as.character(x[1, "feature"]), "composite"),
id = ifelse(length(unique(x[, "id"])) == 1, as.character(x[1, "id"]), cid(i)),
start = min(x[, "start"]),
end = max(x[, "end"]),
group = ifelse(length(unique(x[, "group"])) == 1, as.character(x[1, "group"]), cid(i)),
seqnames = as.vector(nrow(x) == elements[x[1, "seqnames"]])
)
}
})
newVals <- rbind(newVals, as.data.frame(do.call(rbind, tmpVals), stringsAsFactors = FALSE))
merged <- as.logical(newVals$seqnames)
grange <- GRanges(
seqnames = chromosome(GdObject), strand = newVals[, "strand"],
ranges = IRanges(start = as.integer(newVals[, "start"]), end = as.integer(newVals[, "end"]))
)
cnMatch <- match(c(colnames(values(GdObject)), "gdensity"), colnames(newVals))
mcols(grange) <-
if (any(is.na(cnMatch))) newVals[, setdiff(colnames(newVals), c("strand", "start", "end", "seqnames"))] else newVals[, cnMatch]
} else {
grange2 <- GRanges(seqnames = chromosome(GdObject), strand = strand(grange), ranges = ranges(grange))
mcols(grange2) <- mcols(grange)
grange <- grange2
}
return(list(range = grange, needsRestacking = needsRestacking, split = annoSplit, merged = merged, offset = length(annoSplit)))
}
## For AnnotationTracks we need to collapse the all regions along with the additional annotation.
## For GeneRegionTracks we essentially need to do the same thing as for AnnotationTracks, however the additional annotation columns
## are quite different. We do this in multiple turn with increasing levels of complexity:
## 1.) merge all individual items within a group that can no longer be separated
## 2.) merge overlapping groups with just a single remaining item (optional, if mergeGroups==TRUE)
#' @describeIn AnnotationTrack-class preprocess the track before plotting.
#' This will collapse overlapping track items based on the available resolution
#' and increase the width and height of all track objects to a minimum value
#' to avoid rendering issues. See collapsing for details.
#' @keywords internal
setMethod(
"collapseTrack", signature(GdObject = "AnnotationTrack"),
function(GdObject, diff = .pxResolution(coord = "x"), xrange) {
## We first add the original unmodified GdObject as a display parameter to be able to reference back if we ever need to
displayPars(GdObject) <- list(".__OriginalGdObject" = .deepCopyPars(GdObject))
collapse <- .dpOrDefault(GdObject, "collapse", TRUE)
min.width <- .dpOrDefault(GdObject, "min.width", 2)
min.distance <- .dpOrDefault(GdObject, "min.distance", 2)
minXDist <- max(0, ceiling(min.distance * diff))
r <- ranges(GdObject)
## Compute native coordinate equivalent to 1 pixel and resize
rNew <- .resize(r, min.width, diff)
needsRestacking <- any(r != rNew)
r <- rNew
## Collapse all items within a group to a single meta-item (if collapseTranscripts is set)
if (is(GdObject, "GeneRegionTrack")) {
ctrans <- .dpOrDefault(GdObject, "collapseTranscripts", FALSE)
if (is.logical(ctrans) && ctrans) {
ctrans <- "gene"
}
switch(ctrans,
"gene" = {
newVals <- unlist(endoapply(split(values(r), paste(gene(GdObject), strand(GdObject))), head, 1))
newVals$exon <- NA
newVals$feature <- "merged"
newVals$transcript <- newVals$gene
r <- unlist(range(split(r, gene(GdObject))))
mcols(r) <- newVals
GdObject@range <- r
},
"longest" = {
r <- unlist(endoapply(split(r, gene(GdObject)), function(x) {
xs <- split(x, x$transcript)
xs[[which.max(vapply(xs, function(y) abs(diff(as.numeric(as.data.frame(ranges(range(y)))[, c("start", "end")]))), FUN.VALUE = numeric(1L)))]]
}))
GdObject@range <- r
},
"shortest" = {
r <- unlist(endoapply(split(r, gene(GdObject)), function(x) {
xs <- split(x, x$transcript)
xs[[which.min(vapply(xs, function(y) abs(diff(as.numeric(as.data.frame(ranges(range(y)))[, c("start", "end")]))), FUN.VALUE = numeric(1L)))]]
}))
GdObject@range <- r
},
"meta" = {
newVals <- unlist(endoapply(split(values(r), paste(gene(GdObject), strand(GdObject))), head, 1))
newVals$feature <- "merged"
newVals$transcript <- newVals$gene
rtmp <- reduce(split(r, paste(gene(GdObject), strand(GdObject))))
newVals <- newVals[rep(seq_along(rtmp), elementNROWS(rtmp)), ]
newVals$exon <- paste("merged_exon_", unlist(lapply(elementNROWS(rtmp), function(x) seq(1, x)), use.names = FALSE), sep = "")
r <- unlist(rtmp)
mcols(r) <- newVals
GdObject@range <- r
}
)
}
## Collapse overlapping ranges (less than minXDist space between them) and process the annotation data
if (collapse) {
## Merge all items in those groups for which no individual items can be separated
elements <- table(group(GdObject))
rr <- GRanges(seqnames = as.character(group(GdObject)), ranges = IRanges(start = start(r), end = end(r)), strand = strand(r))
mcols(rr) <- mcols(r)
rr <- sort(unique(rr))
mergedAnn <- .collapseAnnotation(rr, minXDist, elements, GdObject)
needsRestacking <- needsRestacking || mergedAnn$needsRestacking
## Now we take a look whether there are any groups that could be merged (if mergeGroups is TRUE)
if (.dpOrDefault(GdObject, "mergeGroups", FALSE) && any(mergedAnn$merged)) {
rr <- sort(mergedAnn$range[mergedAnn$merged])
strand(rr) <- "*"
mergedAnn2 <- .collapseAnnotation(rr, minXDist, elements, GdObject, mergedAnn$offset)
needsRestacking <- needsRestacking || mergedAnn2$needsRestacking
mergedAnn$range <- c(mergedAnn$range[!mergedAnn$merged], mergedAnn2$range)
}
r <- mergedAnn$range
}
## Reconstuct the track object and return
GdObject@range <- r
## if(needsRestacking)
GdObject <- setStacks(GdObject)
return(GdObject)
}
)
## Subset --------------------------------------------------------------------
## In order to keep the grouping information for track regions in the clipped areas we have to
## keep all group elements that overlap with the range. We still want to record the requested
## ranges in the internal '.__plottingRange' display parameter.
#' @describeIn AnnotationTrack-class subset a `AnnotationTrack` by coordinates
#' and sort if necessary.
#' @export
setMethod("subset", signature(x = "AnnotationTrack"), function(x, from = NULL, to = NULL, sort = FALSE, stacks = FALSE, use.defaults = TRUE, ...) {
## Subset to a single chromosome first
lx <- length(x)
csel <- seqnames(x) != chromosome(x)
if (any(csel)) {
x <- x[!csel]
}
if (length(x)) {
## Nothing to do if everything is within the range
granges <- unlist(range(split(ranges(x), group(x))))
ranges <- if (use.defaults) {
.defaultRange(x, from = from, to = to)
} else {
c(
from = ifelse(is.null(from), min(start(granges)) - 1, from),
to = ifelse(is.null(to), max(end(granges)) + 1, to)
)
}
if (!(any(end(x) < ranges["from"] | start(x) > ranges["to"]))) {
if (stacks) {
x <- setStacks(x)
}
return(x)
}
## Now remove everything except for the overlapping groups by first subselecting all groups in the range...
gsel <- names(granges)[subjectHits(findOverlaps(GRanges(seqnames = chromosome(x), ranges = IRanges(min(ranges), max(ranges))), granges))]
x <- x[group(x) %in% gsel]
if (sort) {
x <- x[order(range(x)), ]
}
if (stacks) {
x <- setStacks(x)
}
displayPars(x) <- list(".__plottingRange" = ranges)
}
if (length(x) != lx) {
x <- .computeGroupRange(x)
}
return(x)
})
## ReferenceDataTracks need to stream the data from file and then pass the results on to the next method
#' @describeIn AnnotationTrack-class subset a `ReferenceAnnotationTrack` by
#' coordinates and sort if necessary.
#' @export
setMethod("subset", signature(x = "ReferenceAnnotationTrack"), function(x, from, to, chromosome, ...) {
## We only need to reach out into the referenced file once if the range is already contained in the object
if (missing(from) || is.null(from) || missing(to) || is.null(to)) {
stop("Need both start and end location to subset a ReferenceAnnotationTrack")
}
if (missing(chromosome) || is.null(chromosome)) {
chromosome <- Gviz::chromosome(x)
}
subRegion <- GRanges(seqnames = chromosome[1], ranges = IRanges(start = from, end = to))
if (length(ranges(x)) == 0 || all(overlapsAny(ranges(x), subRegion))) {
cMap <- .resolveColMapping(x@stream(x@reference, subRegion), x@args, x@mapping)
x@range <- .buildRange(cMap$data, args = cMap$args, defaults = x@defaults, trackType = "AnnotationTrack")
chromosome(x) <- chromosome[1]
}
return(callNextMethod(x = x, from = from, to = to, drop = FALSE, ...))
})
## Position ------------------------------------------------------------------
## DrawGrid ------------------------------------------------------------------
setMethod("drawGrid", signature(GdObject = "AnnotationTrack"), function(GdObject, from, to) {
if (.dpOrDefault(GdObject, "grid", FALSE)) {
pushViewport(dataViewport(xData = c(from, to), extension = c(0, 0), yData = 0:1, clip = TRUE))
panel.grid(
h = 0, v = .dpOrDefault(GdObject, "v", -1),
col = .dpOrDefault(GdObject, "col.grid", "#e6e6e6"), lty = .dpOrDefault(GdObject, "lty.grid", 1),
lwd = .dpOrDefault(GdObject, "lwd.grid", 1)
)
popViewport(1)
}
})
## DrawAxis ------------------------------------------------------------------
## DrawGD --------------------------------------------------------------------
## Draw gene models as found in AnnotationTracks or GeneRegionTracks
##
## Calculate all coordinates and values for the individual stacks first, and append to the
## input list 'toDraw'. This allows us to plot the whole track at once, making use of grid's
## vectorization. Without this tweak, every stack would have to be drawn individually, which
## can be painfully slow...
## Because all of the y coordinates are calculated in the coordinate system of the current
## viewport for the respective stack we have to add offset values.
## Compute the coordinates and colors for all track items (e.g. exons in a GeneRegionTrack)
.boxes <- function(GdObject, offsets) {
ylim <- c(0, 1)
h <- diff(ylim)
middle <- mean(ylim)
sh <- max(0, min(h, .dpOrDefault(GdObject, "stackHeight", 0.75)))
space <- (h - (h * sh)) / 2
shift <- switch(.dpOrDefault(GdObject, "stackJust", "middle"),
"top" = space,
"bottom" = -space,
0
)
if (inherits(GdObject, "GeneRegionTrack")) {
thinBox <- .dpOrDefault(GdObject, "thinBoxFeature", .THIN_BOX_FEATURES)
space <- ifelse(feature(GdObject) %in% thinBox, space + ((middle -
space) / 2), space)
}
shape <- .dpOrDefault(GdObject, "shape", "arrow")
color <- .getBiotypeColor(GdObject)
id <- identifier(GdObject, type = .dpOrDefault(
GdObject, ifelse(is(GdObject, "GeneRegionTrack"), "exonAnnotation", "featureAnnotation"),
"lowest"
))
sel <- grepl("\\[Cluster_[0-9]*\\]", id)
id[sel] <- sprintf(
"%i merged\n%s", as.integer(.getAnn(GdObject, "density")[sel]),
ifelse(class(GdObject) %in% c("AnnotationTrack", "DetailsAnnotationTrack"), "features", "exons")
)
yy1 <- ylim[1] + space + offsets + shift
yy2 <- ylim[2] - space + offsets + shift
boxes <- data.frame(
cx1 = start(GdObject), cy1 = yy1, cx2 = start(GdObject) + width(GdObject), cy2 = yy2,
fill = color, strand = strand(GdObject), text = id, textX = start(GdObject) + (width(GdObject) / 2), textY = middle + offsets,
.getImageMap(cbind(start(GdObject), yy1, end(GdObject), yy2)),
start = start(GdObject), end = end(GdObject), values(GdObject), exonId = id, origExonId = identifier(GdObject, type = "lowest"),
stringsAsFactors = FALSE
)
rownames(boxes) <- if (.transcriptsAreCollapsed(GdObject)) {
sprintf("uid%i", seq_along(identifier(GdObject)))
} else {
make.unique(identifier(GdObject, type = "lowest"))
}
return(boxes)
}
## Compute the coordinates for the bars connecting grouped items and the group labels
#' @importFrom grDevices col2rgb hsv rgb2hsv
.barsAndLabels <- function(GdObject) {
bins <- stacks(GdObject)
stacks <- max(bins)
res <- .pxResolution(coord = "x")
gp <- group(GdObject)
if (!is.factor(gp)) {
gp <- factor(gp, levels=unique(gp))
}
grpSplit <- split(range(GdObject), gp)
grpRanges <- unlist(range(grpSplit))
needBar <- vapply(grpSplit, length, FUN.VALUE = numeric(1L)) > 1 & width(grpRanges) > res
## If we draw the bar from start to end of the range we sometimes see little overlaps that extend beyond the first or last item.
## In order to fix this, we just subtract the equivalent of min.width pixels from both ends of each group range
min.swidth <- res * .dpOrDefault(GdObject, "min.width", 2)
nstart <- start(grpRanges[needBar]) + min.swidth
nend <- end(grpRanges[needBar]) - min.swidth
sel <- (nend - nstart) > 0
start(grpRanges[needBar][sel]) <- nstart[sel]
end(grpRanges[needBar][sel]) <- nend[sel]
strand <- vapply(split(strand(GdObject), gp), function(x) {
tmp <- unique(x)
if (length(tmp) > 1) "*" else tmp
}, FUN.VALUE = character(1L))
yloc <- vapply(split((stacks - bins) + 1, gp), function(x) unique(x), FUN.VALUE = numeric(1L)) + 0.5
color <- if (length(grep("__desatCol", values(GdObject)$feature[1]))) {
.dpOrDefault(GdObject, "fill", .DEFAULT_FILL_COL)
} else {
vapply(split(.getBiotypeColor(GdObject), gp), head, 1, FUN.VALUE = character(1L))
}
bars <- data.frame(
sx1 = start(grpRanges)[needBar], sx2 = end(grpRanges)[needBar], y = yloc[needBar], strand = strand[needBar],
col = color[needBar], stringsAsFactors = FALSE
)
labs <- .dpOrDefault(GdObject, ".__groupLabels")[names(grpRanges)]
if (!is.null(labs)) {
lsel <- grepl("\\[Cluster_[0-9]*\\]", labs)
if (any(lsel)) {
gdens <- as.integer(vapply(split(.getAnn(GdObject, "gdensity"), gp), head, 1, FUN.VALUE = character(1L)))
labs[lsel] <- sprintf(
"%i merged %s ", gdens[lsel],
ifelse(class(GdObject) %in% c("AnnotationTrack", "DetailsAnnotationTrack"), "groups", "transcript models")
)
}
just <- .dpOrDefault(GdObject, "just.group", "left")
rev <- .dpOrDefault(GdObject, "reverseStrand", FALSE)
sizes <- .dpOrDefault(GdObject, ".__groupLabelWidths")[names(grpRanges), , drop = FALSE]
pr <- .dpOrDefault(GdObject, ".__plottingRange", data.frame(from = min(start(GdObject)), to = max(end(GdObject))))
## grpRangesCut <- restrict(grpRanges, start=as.integer(pr["from"]), end=as.integer(pr["to"]))
switch(just,
"left" = {
cx <- if (!rev) start(grpRanges) - sizes$after else end(grpRanges) + sizes$after
cy <- yloc
algn <- c("right", "center")
},
"right" = {
cx <- if (!rev) end(grpRanges) + sizes$after else start(grpRanges) - sizes$after
cy <- yloc
algn <- c("left", "center")
},
"above" = {
## cx <- start(grpRangesCut) + width(grpRangesCut)/2
## indx <- which(seq_len(length(grpRanges)) %in% queryHits(findOverlaps(grpRanges, grpRangesCut)))
## cx <- cx[indx] # needs revision
## cy <- yloc[indx] + 0.5
## labs <- labs[indx]
cx <- start(grpRanges) + width(grpRanges) / 2
cy <- yloc + 0.5
algn <- c("center", "top")
},
"below" = {
## cx <- start(grpRangesCut) + width(grpRangesCut)/2
## indx <- which(seq_len(length(grpRanges)) %in% queryHits(findOverlaps(grpRanges, grpRangesCut)))
## cx <- cx[indx] # needs revision
## cy <- yloc[indx] - 0.5
## labs <- labs[indx]
cx <- start(grpRanges) + width(grpRanges) / 2
cy <- yloc - 0.5
algn <- c("center", "bottom")
},
stop(sprintf("Unknown label justification '%s'", just))
)
labels <- data.frame(txt = labs, x = cx, y = cy, stringsAsFactors = FALSE)
} else {
labels <- algn <- NA
}
return(list(bars = bars, labels = labels, align = algn))
}
#' @describeIn AnnotationTrack-class plot the object to a graphics device.
#' The return value of this method is the input object, potentially updated
#' during the plotting operation. Internally, there are two modes in which the
#' method can be called. Either in 'prepare' mode, in which case no plotting is
#' done but the object is preprocessed based on the available space, or in
#' 'plotting' mode, in which case the actual graphical output is created.
#' Since subsetting of the object can be potentially costly, this can be
#' switched off in case subsetting has already been performed before or
#' is not necessary.
#'
#' @export
setMethod("drawGD", signature("AnnotationTrack"), function(GdObject, minBase, maxBase, prepare = FALSE, subset = TRUE, ...) {
debug <- .dpOrDefault(GdObject, "debug", FALSE)
if ((is.logical(debug) && debug) || debug == "prepare") {
browser()
}
imageMap(GdObject) <- NULL
if (!length(GdObject)) {
return(invisible(GdObject))
}
## In prepare mode we need to make sure that the stacking information is updated from the optional display parameter (by calling
## the StackedTrack drawGD method) and also perform the collapsing of track items which could potentially lead to re-stacking.
if (prepare) {
GdObject <- callNextMethod(GdObject, ...)
bins <- stacks(GdObject)
stacks <- max(bins)
## We need to collapse the track object based on the current screen resolution (note that this may trigger re-stacking)
pushViewport(dataViewport(xData = c(minBase, maxBase), extension = 0, yscale = c(1, stacks + 1), clip = TRUE))
GdObject <- collapseTrack(GdObject, diff = .pxResolution(coord = "x"), xrange = c(minBase, maxBase))
popViewport(1)
return(invisible(GdObject))
}
if ((is.logical(debug) && debug) || debug == "draw") {
browser()
}
## If there are too many stacks for the available device resolution we cast an error, otherwise we set the viewport
bins <- stacks(GdObject)
stacks <- max(bins)
rev <- .dpOrDefault(GdObject, "reverseStrand", FALSE)
xscale <- if (!rev) c(minBase, maxBase) else c(maxBase, minBase)
yscale <- if (!.dpOrDefault(GdObject, "reverseStacking", FALSE)) c(1, stacks + 1) else c(stacks + 1, 1)
pushViewport(dataViewport(xscale = xscale, extension = 0, yscale = yscale, clip = TRUE))
res <- .pxResolution(coord = "x")
curVp <- vpLocation()
if (curVp$size["height"] / stacks < .dpOrDefault(GdObject, "min.height", 3)) {
stop("Too many stacks to draw. Either increase the device size or limit the drawing to a smaller region.")
}
## We adjust the color saturation to indicate overplotting if necessary
if (.dpOrDefault(GdObject, "showOverplotting", FALSE)) {
dens <- as.numeric(values(GdObject)$density)
if (length(unique(dens)) != 1) {
minSat <- max(0.25, 1 / max(dens))
minDens <- min(dens)
rDens <- diff(range(dens))
saturation <- minSat + ((dens - minDens) / rDens / (1 / (1 - minSat)))
bc <- unique(.getBiotypeColor(GdObject))
baseCol <- rgb2hsv(col2rgb(bc))
desatCols <- unlist(lapply(saturation, function(x) hsv(baseCol[1, ], x, baseCol[3, ])))
names(desatCols) <- paste(unique(feature(GdObject)), rep(dens, each = length(bc)), sep = "_")
feature(GdObject) <- paste(feature(GdObject), dens, sep = "_")
desatCols <- desatCols[unique(names(desatCols))]
displayPars(GdObject) <- as.list(desatCols)
}
}
## Now we can pre-compute all the coordinates and settings for the elements to be drawn, ...
box <- .boxes(GdObject, (stacks - bins) + 1)
barsAndLab <- .barsAndLabels(GdObject)
bar <- barsAndLab$bars
bartext <- barsAndLab$labels
## ... get all the necessary display parameters
shape <- .dpOrDefault(GdObject, "shape", "arrow")
col.line <- .dpOrDefault(GdObject, "col.line")[1]
border <- .dpOrDefault(GdObject, "col")[1]
if (is.null(border)) {
border <- ifelse(is(GdObject, "GeneRegionTrack"), NA, "transparent")
}
lwd <- .dpOrDefault(GdObject, "lwd", 2)
lty <- .dpOrDefault(GdObject, "lty", 1)
alpha <- .dpOrDefault(GdObject, "alpha", 1)
rotation <- .dpOrDefaultFont(GdObject, "rotation", "item", 0)
rotation.group <- .dpOrDefaultFont(GdObject, "rotation", "group", 0)
just <- .dpOrDefault(GdObject, "just.group", "left")
## ... and finally draw whatever is needed
if (nrow(box) > 0) {
## If we want to place a label on top or below the ranges we need to know how much size that will take up
## and adjust the plotting shapes accordingly
drawLabel <- .dpOrDefault(GdObject, ".__hasAnno", FALSE) && !is.null(bartext) && !anyNA(bartext) &&
nrow(bartext) > 0 && stacking(GdObject) != "dense"
bs <- as.vector((stacks - bins) + 1)
if (drawLabel && just %in% c("above", "below")) {
labelHeights <- max(.getStringDims(GdObject, bartext$txt, subtype = "group")$height)
labelSpace <- as.numeric(convertHeight(unit(2, "points"), "native"))
avSpace <- (1 - max(box$cy2 - box$cy1)) / 2
vadjust <- max(0, (labelHeights + (2 * labelSpace)) - avSpace)
bh <- 1 - (avSpace * 2)
bhNew <- bh - vadjust
sfac <- bhNew / bh
if (just == "above") {
bar$y <- bar$y - (vadjust / 2)
box$textY <- box$textY - (vadjust / 2)
bartext$y <- bartext$y - labelSpace
box$cy1 <- bs + ((box$cy1 %% 1 - avSpace) * sfac) + avSpace
box$cy2 <- bs + ((box$cy2 %% 1 - avSpace) * sfac) + avSpace
} else {
bar$y <- bar$y + (vadjust / 2)
box$textY <- box$textY + (vadjust / 2)
bartext$y <- bartext$y + labelSpace
box$cy1 <- bs + ((box$cy1 %% 1 - avSpace) * sfac) + avSpace + vadjust
box$cy2 <- bs + ((box$cy2 %% 1 - avSpace) * sfac) + avSpace + vadjust
}
}
## Plotting of the (arrow)bar
if (nrow(bar) > 0) {
.arrowBar(bar$sx1, bar$sx2,
y = bar$y, bar$strand, box[, seq_len(4), drop = FALSE],
W = .dpOrDefault(GdObject, "arrowFeatherWidth", 3), D = .dpOrDefault(GdObject, "arrowFeatherDistance", 20),
col = if (is.null(col.line)) bar$col else rep(col.line, length(bar$col)), lwd = lwd, lty = lty,
alpha = alpha, barOnly = (!"smallArrow" %in% .dpOrDefault(GdObject, "shape", "box") || stacking(GdObject) == "dense"),
diff = res, min.height = .dpOrDefault(GdObject, "min.height", 3)
)
}
## Plotting of the boxes
box$col <- if (is.na(border)) box$fill else border
if ("box" %in% shape || ("smallArrow" %in% shape && !("arrow" %in% shape || "ellipse" %in% shape))) {
.filledBoxes(box, lwd = lwd, lty = lty, alpha = alpha)
}
## Plotting of the ellipses
if ("ellipse" %in% shape) {
ellCoords <- .box2Ellipse(box)
grid.polygon(
x = ellCoords$x1, y = ellCoords$y1, id = ellCoords$id,
gp = gpar(col = if (is.na(border)) box$fill else border, fill = box$fill, lwd = lwd, lty = lty, alpha = alpha),
default.units = "native"
)
}
## Plotting of the filled arrows
if ("arrow" %in% shape && !"box" %in% shape) {
.filledArrow(box, lwd = lwd, lty = lty, alpha = alpha, min.width = 4 * res, max.width = .dpOrDefault(GdObject, "arrowHeadMaxWidth", 40) * res)
}
## Plotting of the filled arrows with fixed head size
if ("fixedArrow" %in% shape && !"box" %in% shape) {
.filledArrow(box,
lwd = lwd, lty = lty, alpha = alpha, min.width = 4 * res, absoluteWidth = TRUE,
W = .dpOrDefault(GdObject, "arrowHeadWidth", 30) * res
)
}
## Plotting of the item labels
if (.dpOrDefault(GdObject, "showFeatureId", FALSE)) {
# grid.text(str2expression(box$text), box$textX, box$textY,
grid.text(box$text, box$textX, box$textY,
rot = rotation, gp = .fontGp(GdObject, subtype = "item"),
default.units = "native", just = c("center", "center")
)
}
## Plotting of the group labels
if (drawLabel) {
grid.text(bartext$txt, bartext$x, bartext$y,
rot = rotation.group, gp = .fontGp(GdObject, subtype = "group"),
default.units = "native", just = barsAndLab$align
)
}
}
popViewport(1)
## Finally we set up the image map
## FIXME: we may want to record the merging information here
im <- if (!is.null(box)) {
coords <- as.matrix(box[, c("x1", "y1", "x2", "y2"), drop = FALSE])
restCols <- setdiff(colnames(box), c("x1", "x2", "y1", "y2", "cx1", "cx2", "cy1", "cy2", "textX", "textY"))
tags <- lapply(restCols, function(x) {
tmp <- as.character(box[, x])
names(tmp) <- rownames(coords)
tmp
})
names(tags) <- restCols
tags$title <- identifier(GdObject)
ImageMap(coords = coords, tags = tags)
} else {
NULL
}
imageMap(GdObject) <- im
return(invisible(GdObject))
})
## DrawGD DetailsAnnotationTrack ---------------------------------------------
## Create a data.frame with the distinct details function arguments (like start, end, ...)
.buildArgsDf <- function(GdObject) {
groupDetails <- .dpOrDefault(GdObject, "groupDetails", FALSE)
rr <- if (groupDetails) unlist(range(split(ranges(GdObject), group(GdObject)))) else ranges(GdObject)
args <- data.frame(
start = as.integer(start(rr)), end = as.integer(end(rr)), strand = as.character(strand(rr)),
chromosome = as.character(seqnames(rr)),
identifier = as.character(if (groupDetails) names(rr) else identifier(GdObject, type = "lowest")),
stringsAsFactors = FALSE
)
return(args)
}
#' @describeIn AnnotationTrack-class plot the object to a graphics device.
#' The return value of this method is the input object, potentially updated
#' during the plotting operation. Internally, there are two modes in which the
#' method can be called. Either in 'prepare' mode, in which case no plotting is
#' done but the object is preprocessed based on the available space, or in
#' 'plotting' mode, in which case the actual graphical output is created.
#' Since subsetting of the object can be potentially costly, this can be
#' switched off in case subsetting has already been performed before or
#' is not necessary.
#'
#' @export
setMethod(
"drawGD", signature("DetailsAnnotationTrack"),
function(GdObject, minBase, maxBase, prepare = FALSE, ...) {
debug <- .dpOrDefault(GdObject, "debug", FALSE)
if ((is.logical(debug) && debug) || debug == "prepare") {
browser()
}
adf <- .buildArgsDf(GdObject)
args <- .dpOrDefault(GdObject, "detailsFunArgs", fromPrototype = TRUE)
groupDetails <- .dpOrDefault(GdObject, "groupDetails", FALSE)
if (prepare) {
GdObject <- callNextMethod()
GdObject <- GdObject[order(start(GdObject))]
indices <- if (groupDetails) seq_len(length(unique(group(GdObject)))) else seq_len(length(GdObject))
xscale <- if (!.dpOrDefault(GdObject, "reverseStrand", FALSE)) c(minBase, maxBase) else c(maxBase, minBase)
pushViewport(viewport(xscale = xscale))
select <- lapply(indices, function(i) {
warn <- FALSE
iargs <- as.list(adf[i, ])
iargs$index <- i
iargs$GdObject <- GdObject
iargs$GdObject.original <- .dpOrDefault(GdObject, ".__OriginalGdObject", GdObject)
args <- c(args[setdiff(names(args), names(iargs))], iargs)
res <- do.call(GdObject@selectFun, args)
if (length(res) != 1 || !is.logical(res) || is.na(res)) {
res <- TRUE
warn <- TRUE
}
c(res = res, warn = warn)
})
select <- do.call(rbind, select)
if (any(select[, "warn"])) {
warning("The result of function 'selectFun' has to be a single logical value. Forcing the value to 'TRUE'")
}
select <- select[, "res"]
popViewport(1)
displayPars(GdObject) <- list(".__select" = select)
return(invisible(GdObject))
}
if ((is.logical(debug) && debug) || debug == "draw") {
browser()
}
n <- length(GdObject)
col <- rep(.dpOrDefault(GdObject, "detailsConnector.col", fromPrototype = TRUE), n)[seq_len(n)]
lty <- rep(.dpOrDefault(GdObject, "detailsConnector.lty", fromPrototype = TRUE), n)[seq_len(n)]
lwd <- rep(.dpOrDefault(GdObject, "detailsConnector.lwd", fromPrototype = TRUE), n)[seq_len(n)]
pch <- rep(.dpOrDefault(GdObject, "detailsConnector.pch", fromPrototype = TRUE), n)[seq_len(n)]
cex <- rep(.dpOrDefault(GdObject, "detailsConnector.cex", fromPrototype = TRUE), n)[seq_len(n)]
border.lty <- rep(.dpOrDefault(GdObject, "detailsBorder.lty", fromPrototype = TRUE), n)[seq_len(n)]
border.lwd <- rep(.dpOrDefault(GdObject, "detailsBorder.lwd", fromPrototype = TRUE), n)[seq_len(n)]
border.col <- rep(.dpOrDefault(GdObject, "detailsBorder.col", fromPrototype = TRUE), n)[seq_len(n)]
border.fill <- rep(.dpOrDefault(GdObject, "detailsBorder.fill", fromPrototype = TRUE), n)[seq_len(n)]
minwidth <- .dpOrDefault(GdObject, "details.minWidth", fromPrototype = TRUE)
size <- .dpOrDefault(GdObject, "details.size", fromPrototype = TRUE)
xyratio <- .dpOrDefault(GdObject, "details.ratio", fromPrototype = TRUE)
if (0 >= size || size > 1) {
warning("details.size must be >0 and <1 - reset to 0.5")
size <- 0.5
}
selection <- .dpOrDefault(GdObject, ".__select", rep(TRUE, length(GdObject)))
len <- sum(selection)
bins <- if (!groupDetails) stacks(GdObject) else unlist(lapply(split(stacks(GdObject), group(GdObject)), unique))
stacks <- max(bins)
if (len > 0) {
if (((maxBase - minBase) / len) / .pxResolution(coord = "x") < minwidth) {
warning("too much detail for available space (plot fewer annotation or increase details.minWidth)!")
popViewport(1)
GdObject <- callNextMethod()
return(GdObject)
}
rr <- if (groupDetails) unlist(range(split(ranges(GdObject), group(GdObject)))) else ranges(GdObject)
xloc1 <- (end(rr) - start(rr)) / 2 + start(rr)
yloc1 <- (stacks - (bins - 0.5) + 1)
xloc2 <- ((1 / len * seq_len(len)) - 1 / len + (1 / len * 0.5))
yloc2 <- rep(1, len)
## draw details plots (via user supplied function 'fun')
pushViewport(viewport(height = size, y = 1 - size, just = c(0.5, 0)))
w <- 1
v <- 0
vpl <- vpLocation()
r <- vpl$size["width"] / len / vpl$size["height"]
if (r > xyratio) {
w <- xyratio / r
v <- ((1 / len) - (1 / len * w)) / 2
}
indices <- if (groupDetails) seq_len(length(unique(group(GdObject)))) else seq_len(length(GdObject))
j <- 1
pres <- list()
hasError <- FALSE
for (i in indices[selection]) {
pushViewport(viewport(width = 1 / len * w, x = ((1 / len * j) - 1 / len) + (v), just = c(0, 0.5)))
grid.rect(gp = gpar(col = border.col[i], lwd = border.lwd[i], lty = border.lty[i], fill = border.fill[i]))
iargs <- as.list(adf[i, ])
iargs$index <- i
iargs$GdObject <- GdObject
iargs$GdObject.original <- .dpOrDefault(GdObject, ".__OriginalGdObject", GdObject)
args <- c(args[setdiff(names(args), names(iargs))], iargs)
pres[[as.character(j)]] <- try(do.call(GdObject@fun, args), silent = TRUE)
if (!is.null(pres) && is(pres[[as.character(j)]], "try-error")) {
hasError <- TRUE
grid.segments(x0 = c(0.1, 0.1), x1 = c(0.9, 0.9), y0 = c(0.9, 0.1), y1 = c(0.1, 0.9), gp = gpar(col = "red", lwd = 3))
}
popViewport(1)
j <- j + 1
}
if (hasError) {
warning("There have been errors in the detail plotting function:\n", paste(pres, collapse = "\n"))
}
popViewport(1)
## plot AnnotationTrack and connectors to details
xscale <- if (!.dpOrDefault(GdObject, "reverseStrand", FALSE)) c(minBase, maxBase) else c(maxBase, minBase)
pushViewport(viewport(
xscale = xscale,
yscale = c(1, stacks + 1), clip = FALSE,
height = 1 - size, y = 0, just = c(.5, 0)
))
GdObject <- callNextMethod()
grid.segments(
x0 = unit(xloc1[selection], "native"), x1 = xloc2, y0 = unit(yloc1[selection], "native"),
y1 = yloc2, gp = gpar(col = col, lwd = lwd, lty = lty, cex = cex)
)
grid.points(x = unit(xloc2, "npc"), y = unit(yloc2, "npc"), gp = gpar(col = col, cex = cex), pch = pch)
grid.points(x = unit(xloc1[selection], "native"), y = unit(yloc1[selection], "native"), gp = gpar(col = col, cex = cex), pch = pch)
popViewport(1)
} else {
GdObject <- callNextMethod()
}
return(invisible(GdObject))
}
)
## SetAs ---------------------------------------------------------------------
#' @importFrom rtracklayer GenomicData
setAs(
"AnnotationTrack", "UCSCData",
function(from, to) {
ranges <- range(from)
dcolor <- as.integer(col2rgb(.dpOrDefault(from, "col")))
line <- new("BasicTrackLine",
name = names(from),
description = names(from),
visibility = stacking(from), color = dcolor, itemRgb = TRUE
)
vals <- values(from)
color <- .getBiotypeColor(from)
strand <- as.character(strand(from))
strand[strand == "*"] <- "+"
new("UCSCData", rtracklayer::GenomicData(ranges,
chrom = chromosome(from),
id = gsub(" ", "_", vals$id),
name = gsub(" ", "_", as.character(vals$id)), itemRgb = color,
strand = strand
),
trackLine = line
)
}
)
setAs("GRanges", "AnnotationTrack", function(from, to) AnnotationTrack(range = from))
setAs("GRangesList", "AnnotationTrack", function(from, to) AnnotationTrack(range = from))
## Show ----------------------------------------------------------------------
## A helper function to plot general information about an AnnotationTrack
.annotationTrackInfo <- function(object) {
msg <- sprintf(
paste("| genome: %s\n| active chromosome: %s\n",
"| annotation features: %s",
sep = ""
),
genome(object),
chromosome(object),
length(object)
)
addfeat <- length(object@range) - length(object)
if (addfeat > 0) {
msg <- c(msg, .addFeatInfo(object, addfeat), "Call chromosome(obj) <- 'chrId' to change the active chromosome")
}
return(paste(msg, collapse = "\n"))
}
## We have to show the name, genome and currently active chromosome, and, if more ranges are available on additional
## chromosomes some information about that
#' @describeIn AnnotationTrack-class Show method.
#' @export
setMethod("show", signature(object = "AnnotationTrack"), function(object) {
cat(sprintf("AnnotationTrack '%s'\n%s\n", names(object), .annotationTrackInfo(object)))
})
#' @describeIn AnnotationTrack-class Show method.
#' @export
setMethod("show", signature(object = "ReferenceAnnotationTrack"), function(object) {
.referenceTrackInfo(object, "ReferenceAnnotationTrack")
})
