#' @exportMethod with
#' @exportMethod within
#' @importFrom data.table data.table
#'
#' @section Slots:
#' \describe{
#' \item{data}{length(.Object) length list containing genomic data (e.g. GRanges, GrangesLists, RleLists, path to ucsc file or ffTrack file on disk)}
#' \item{seqinfo}{Seqinfo object}
#' \item{colormap}{length(.Object) length named list of named vectors whose entry i specifies the colormap for the meta data field of object entry i (specified by the name) and maps unique value of that data field to colors (specified by the named vector)}
#' \item{edges}{list of data.frames of length length(.Object) which has columns $from, $to, and optional fields $col, $lwd, and $lty to specify splined edges joining data items in the corresponding track}
#' \item{formatting}{\code{data.frame} holding all of the formatting options}
#' \item{mdata}{\code{matrix} holding interaction data between genomic loci (e.g. hiC)}
#' \item{vars}{List of variants associated with track}
#' }
#'
#' @name gTrack-class
#' @rdname gTrack-class
#' @description
#' S4 class for \code{gTrack}
#'
#' Class \code{gTrack} defines a subsettable object that wraps formatting information around
#' genomic data (GRanges, GRangesList, RleLists, UCSC formats, ffTrack) and can be displayed
#' in a "genome browser" style plot.
#'
#' @exportClass gTrack
#' @author Marcin Imielinski
#' @importFrom methods setClass setGeneric setMethod setRefClass
#' @import rtracklayer
#' @importFrom gUtils grl.unlist si2gr grbind gr.string gr.fix grl.pivot gr.findoverlaps gr.flatten gr.chr gr.match
setClass('gTrack', representation(data = 'list', mdata= 'list', seqinfo = 'Seqinfo', formatting = 'data.frame', colormap = 'list', edges = 'list', vars = 'list'))
#setClass('trackData', contains = "gTrack") ## for legacy, backwards compatibility with old trackData class
setMethod('initialize', 'gTrack', function(.Object, data, mdata, edges, vars, colormaps, seqinfo, y.field, yaxis, format, ...) ## only place NON formatting fields here. The rest passed with ...
{
.Object@data <- listify(data, GRanges)
if (any(ix <- !(sapply(.Object@data, class) %in% c('GRanges', 'GRangesList', 'CompressedGRangesList', 'character', 'RleList', 'ffTrack'))))
stop('check input: gTrack objects can only be defined around GRanges, GRangesLists, RleLists, ffTrack, file paths to .rds files of the latter object types, or file paths to UCSC format files')
## replicate data and/or y.field to match up data and y.field lengths
if (length(y.field) != length(.Object@data))
{
repdt = data.table(data.id = 1:length(.Object@data), y.field = y.field)
.Object@data = .Object@data[repdt$data.id]
y.field = repdt$y.field
}
.Object@mdata <- listify(mdata, matrix, length(.Object@data))
for (i in 1:length(.Object))
if (!is.matrix(.Object@mdata[[i]]) & !inherits(.Object@mdata[[i]], 'Matrix'))
stop('optional arg mdata be either empty matrix or a square matrix of same dimensions as data GRanges')
else if (!identical(.Object@mdata[[i]], matrix()))
{
if (!identical(dim(.Object@mdata[[i]]), rep(length(.Object@data[[i]]), 2)))
stop('mdata for each entry must be a square matrix of same dimensions as data GRanges')
if (inherits(.Object@mdata[[i]], 'Matrix'))
.Object@mdata[[i]] = (Matrix::t(.Object@mdata[[i]]) + .Object@mdata[[i]])/2
else
.Object@mdata[[i]] = (t(.Object@mdata[[i]]) + .Object@mdata[[i]])/2
}
# .Object@edges <- listify(edges, data.frame, length(.Object@data)) ## listify not working here i.e. for concatenating objects with @edges field
if (is.null(edges))
.Object@edges = rep(list(data.frame()), length(.Object@data))
else if (!is.list(edges) | inherits(edges, 'data.frame'))
.Object@edges = list(edges)
else
.Object@edges = edges
.Object@vars <- listify(vars, list, length(.Object@data))
null.ix = sapply(.Object@data, is.null)
if (any(null.ix))
.Object@data[null.ix] = list(GRanges())
edges.df = sapply(.Object@edges, is.data.frame)
edges.mat = sapply(.Object@edges, function(x) is.array(x) | inherits(x, 'Matrix'))
if (any(edges.df))
.Object@edges[edges.df] = lapply(.Object@edges[edges.df], function(x)
{
if (nrow(x)>0)
{
if (!all(c('from', 'to') %in% names(x)))
stop('edges data frame missing $to and $from columns')
if (data.table::is.data.table(x))
x = as.data.frame(x)
x
}
else
data.frame()
})
if (any(edges.mat))
.Object@edges[edges.mat] = lapply(which(edges.mat), function(x)
{
y = .Object@edges[[x]]
if (nrow(y) != length(.Object@data[[x]]) & ncol(y) != length(.Object@data[[x]]))
warning('dimensions of matrix do not match length of track data item. Adjacency matrix of edges should be square and have as many rows as there are input ranges in the corresponding gTrack item')
tmp.edges = Matrix::which(y>0, arr.ind = TRUE)
return(data.frame(from = tmp.edges[,1], to = tmp.edges[,2], lwd = y[tmp.edges]))
})
.Object@colormap <- listify(colormaps, list)
if (length(.Object@colormap)==1)
.Object@colormap = rep(.Object@colormap, length(.Object@data))
if (is.null(names(.Object@colormap)))
names(.Object@colormap) = NA
## convert options to formatting string
if (!is.data.frame(format) && is.na(format))
.Object@formatting = as.data.frame(c(list(y.field = y.field), list(...)), stringsAsFactors=FALSE)
else if (is.data.frame(format) && nrow(format) == length(.Object))
.Object@formatting <- cbind(data.frame(y.field = y.field), format) #formatting(.Object) <- format
else
stop("Expecting formatting (as supplied directly to constructor) to be data.frame of nrow = length(gTrack)")
## figure out seqinfo from provided tracks, make sure that all the seqinfos are compatible
## if only Rle's provided
## JEREMIAH 3/20/16 t.name does not need default
.Object <- get_seqinfo(.Object, seqinfo)
## if single data input but multiple yfields, replicate
if (length(y.field) > 1 && length(.Object@data) == 1) {
new_object = .Object
new_object$y.field = y.field[1]
new_object$yaxis = yaxis[1]
for (i in 2:length(y.field)) {
o = .Object
o$y.field = y.field[i]
o$yaxis = yaxis[i]
new_object = c(new_object, o)
}
.Object = new_object
} else if (!all(is.na(y.field))) {
formatting(.Object)$y.field <- y.field
formatting(.Object)$yaxis <- yaxis
}
else {
formatting(.Object)$y.field <- NA
formatting(.Object)$yaxis <- TRUE
}
if (any(!is.na(y.field))) {
ix <- as.logical(pmax(nchar(formatting(.Object)$name) == 0 , is.na(formatting(.Object)$name),na.rm = TRUE))
formatting(.Object)$name[ix] <- y.field[ix]
}
return(.Object)
})
#' Construct a new \code{gTrack}
#'
#' Arguments described as "formatting" are vectors. They are replicated (if necessary) to match the length of the object.
#' @param data Instance of one of the following objects: (1) GRanges (2) GRangesList (3) ffTrack or
#' (4) character representing Bed, Wig, BigWig, or .rds GRanges file. It can also be a list of the above (if specifying multiplying
#' tracks.
#' @param mdata Square matrix representing values at connections between genomic elements (e.g. chr8:1,000-2,000 to chr3:5,000-6,000). The number of
#' rows (and thus columns) must be equal to the length of the \code{data} element. The value of the (i,j) element of this matrix represents a value
#' for the intersection of \code{data[i]} and \code{data[j]}. Note that all connections are assumed to be symmetric, and only the upper triangle
#' will be read. Displays a \code{warning} if the matrix is non-symmetric
#' @param y.field vector or scalar character referencing meta data field of GRanges or GRangesList to use for "y axis" coordinate when plotting
#' numeric tracks, (note: for a RleList or ffTrack this field is automatically set to "score"), default is NA for non-numeric tracks
#' @param name vector or scalar character specifying name of this track, which will be displayed on label to the left of the track
#' @param height vector or scalar numeric specifying height of track(s) (in relative units)
#' @param gr.labelfield vector or scalar character specifying which GRanges meta data field to use for GRanges label (default "label") (formatting)
#' @param grl.labelfield vector or scalar character specifying which GRanges meta data field to use for GRangesList label (default "label") (formatting)
#' @param grl.colorfield vector or scalar character specifying which GRanges meta data field to use for GRangesList labels color (default "label.color") (formatting)
#' @param grl.cexfield vector or scalar character specifying which GRanges meta data field to use for GRangesList labels cex (default "cex.label") (formatting)
#' @param legend.maxitems Scalar positive integer specifying what is the maximum number of items to include in legend [Inf]
#' @param label.suppress vector or scalar logical flag specifying whether to suppress all GRanges / GRangesList label drawing (formatting)
#' @param label.suppress.gr vector or scalar logical flag specifying whether to suppress GRanges label drawing (formatting)
#' @param label.suppress.grl vector or scalar logical flag specifying whether to suppress GRangesList label drawing (formatting)
#' @param ygap vector or scalar numeric specifying gap between tracks
#' @param stack.gap vector or scalar.numeric specifying x gap between stacking non-numeric GRanges or GrangesLists items in track(s)
#' @param cex.label size of top level label (e.g. GRAngesList item)
#' @param cex.label.gr size of bottom level label (i.e. GRanges in a GrangesListitem, e.g. exons in a gene)
#' @param gr.srt.label rotation of GRanges label
#' @param ywid vector or scalar numeric specifying the y-extent of individual ranges (in local plot coordinates)
#' @param col vector or scalar character specifying static color for track(s), if NA then color is specified by colormaps() property or gr.colorfield or col meta data field of GRanges / GRangesList data object
#' @param border vector or scalar character specifying static border for polygons in track(s), if NA then $border is determine duing gr.colorfield / colormap or meta field $border of GRanges / GRangesList
#' @param max.ranges vector or scalar numeric specifying what is the max number of ranges to plot in a window (formatting)
#' @param angle vector of scalar numeric specifying angle of polygons that represent signed range'
#' items (only relevant if used within chainedTracks object)
#' @param split vector or scalar logical flag specifying whether to split when lifting (only relevant if used wihtin chainedTracks object)
#' @param colormaps length(.Object) length named list of named vectors whose entry i maps uniques value of a data field to colors. The data.field is specified by the name of the list entry, and the unique values / colors are specified by the named vector.
#' @param edges Data frame of columns $from, $to, and optional fields $col, $lwd, and $lty, specifying edges linking data items.
#' Also can be a list of the above if specifying multiple tracks (and must be compatible in length with data arg)
#' @param xaxis.suffix vector or scalar numeric specifying the suffix that will be used in describin x axis coordinates (TODO: move to display) (formatting)
#' @param xaxis.unit vector or scalar numeric specifying the unit that will be used in describing x axis coordinates (TODO: move to display) (formatting)
#' @param xaxis.round vector or scalar non-neg integer specifying number of decimals to round xaxis coordinate labels (formatting)
#' @param xaxis.nticks vector or scalar positive integer specifying how many xaxis ticks to optimally draw (formatting)
#' @param xaxis.label.angle vector or scalar numeric between 0 and 360 specifying angle with which to draw xaxis coordinate labels (formatting)
#' @param xaxis.newline vector or scalar logical specifying whether to draw a newline in the xaxis coordinate labels (formatting)
#' @param xaxis.width Logical scalar specifying whether to add window width to xaxis window labels [TRUE]
#' @param lwd.border vector or scalar integer specifying the thickness of the polygon borders (formatting)
#' @param cex.label vector or scalar numeric specifying the expansion factor of the range labels (formatting)
#' @param hadj.label vector or scalar numeric specifying the horizontal adjustment of the range labels (formatting)
#' @param vadj.label vector or scalar numeric specifying the vertical adjustment of the range labels (formatting)
#' @param ypad vector or scalar numeric between 0 and 1 specifying how much whitespace padding to add within panel (formatting)
#' @param circles vector or scalar logical specifying whether to scatter plot range data (formatting)
#' @param lines vector or scalar logical specifying whether to line plot range data (formatting)
#' @param bars vector or scalar logical specifying whether to bar plot range data (formatting)
#' @param y0.bar vector or scalar numeric specifying where to draw the lower boundary of a bar in a bar plot (only applicable if bars == T) (formatting)
#' @param source.file.chrsub vector or scalar logical specifying whether or not sub "chr" out of any external files (e.g. UCSC style files) (formatting)
#' @param y.grid.col vector or scalar character specifying color of "gridlines" used to specify numeric track data (formatting)
#' @param y.grid.cex vector or scalar non-neg numeric specifying character expansion for y tick / y grid labels (formatting)
#' @param y.grid.lty vector or scalar positive integer specifying line style of y grid lines for numeric tracks (formatting)
#' @param y.grid.lwd vector or scalar positive integer specifying thickness of y grid lines for numeric tracks (formatting)
#' @param y.grid.labx vector or scalar positive integer specifying fraction of xlim left of plot to place y axis labels (formatting)
#' @param yaxis vector or scalar logical specifying whether to print yaxis (formatting)
#' @param yaxis.pretty vector or scalar positive integer specifying how many ticks to optimally draw on yaxis (formatting)
#' @param draw.var vector or scalar logical specifying whether to draw.var for GRanges / GRangesList specifying reads
#' (GRanges must contain $cigar +/- $MD field) (formatting)
#' @param draw.paths vector or scalar logical specifying whether to interpret GRangesLists as "paths" and connect them with a
#' a set of spline curves. (formatting)
#' @param path.col vector or scalar character specifying color of path (only applicable for tracks in which draw.paths = T) (formatting)
#' @param path.col.arrow vector or scalar character specifying color of arrow of path (only applicable for tracks in which draw.paths = T) (formatting)
#' @param path.cex.arrow vector or scalar numeric > 0 specifying expansion factor of arrow of path (only applicable for tracks
#' in which draw.paths = T) (formatting)
#' @param path.stack.y.gap vector or scalar numeric > 0 specifying y stack gap of paths (only applicable for tracks in which draw.paths = T) (formatting)
#' @param path.stack.x.gap vector or scalar numeric > 0 specifying x stack gap for paths (only applicable for tracks in which draw.paths = T) (formatting)
#' @param path.cex.v vector or scalar numeric > 0 specifying vertical bulge of sline in paths (only applicable for tracks in which draw.paths = T) (formatting)
#' @param path.cex.h vector or scalar numeric > 0 specifying horizontal bulge of spline in paths (formatting)
#' (only applicable for tracks in which draw.paths = T) (formatting)
#' @param draw.backbone vector or scalar logical specifying whether to draw "backbone" connecting different items in a GRangesList item (formatting)
#' @param xaxis.cex.tick Scalar numeric specifying expansion factor for axis tick labels (formatting)
#' @param xaxis.ticklen Scalar numeric specifying lengths for axis ticks (formatting)
#' @param gr.cex.label vector or scalar numeric specifying GRanges label character expansion (default is cex.label) (formatting)
#' @param gr.srt.label vector or scalar numeric between 0 and 180 specifying rotation of GRanges labels (formatting)
#' @param sep.lty Scalar integer specifying line style for window separators [2] (dashed line)
#' @param sep.lwd Scalar numeric specifying line thickness for window separators [1]
#' @param sep.bg.col Color (supplied by name or hex code) for the background of the windows [gray95]
#' @param sep.draw Logical allowing separators to be turned off [FALSE]
#' @param cmap.min minimum saturating data value of color map for triangle plot
#' @param cmap.max maximum saturating data value of color map for triangle plot
#' @param labels.suppress whether to suppress labels (both GRangesList and GRanges)
#' @param labels.suppress.gr whether to suppress GRanges labels
#' @param labels.suppress.gr whether to suppress GRangesList labels
#' @param ... additional formatting arguments to gTrack
#'
#' @name gTrack
#' @rdname gTrack-class
#' @export
#' @author Marcin Imielinski
gTrack = function(data = NULL, ##
y.field = NA, ## character field specifying what numeric values field of the underlying gr or grl should be used as y coordinate (otherwise ranges are stacked)
mdata = NULL, ## this is matrix of values for triangle plot
edges = NULL, ## this is a list of data.frames of length numtracks or a data.frame, data.frame has fields $from, $to, and optional fields $col, $lwd, $lty to specify edges joining ranges and their display, can also be an adjacency matrix of dim n x n where n = length(GRanges) or length(GRangesList), the values of the matrix will specify $lwd of the resulting lines
vars = NULL, ## GRangesList of same length as data (if length is 1 ) or list of GRangesList representing variants (output of draw.var)
colormaps = NULL, # (named) list same length as data
height = 10,
ygap = 2,
ylab = '',
stack.gap = 0,
cex.label = 1,
gr.cex.label = cex.label *0.8,
gr.srt.label = 0,
col = NA,
border = NA,
angle = 15,
name = "",
gr.colorfield = NA,
y.quantile = 0.01, ## if y0 or y1 is not specified then will draw between y.quantile and 1-y.quantile of data
y.cap = T, ## whether to cap values at y0, y1 (only relevant if y.field specified)
lwd.border = 1,
hadj.label = 1,
vadj.label = 0.5,
smooth = NA, ## smooth with running mean with this window
round = NA, ## round the output of running mean to this precision
ywid = NA,
ypad = 0,
seqinfo = NA,
circles = FALSE,
lines = FALSE,
bars = FALSE,
draw.paths = FALSE,
path.col = 'black',
path.lwd = 1,
draw.var = FALSE,
triangle = !is.null(mdata),
max.ranges = 5e4, ## parameter to limit max number of ranges to draw on canvas, will downsample to this amount
source.file.chrsub = T, ## if source file has chr for seqnames this will sub it out
y0.bar = NA,
yaxis = !is.na(y.field), # logical whether to print yaxis
yaxis.pretty = 5, # how many ticks to optimally draw on yaxis
yaxis.cex = 1, ## size of text at yaxis
chr.sub = TRUE, ## remove 'chr' from slen if drawing from UCSC style format
edgevars = NA,
gr.labelfield = NA,
grl.labelfield = NA,
grl.colorfield = NA,
grl.cexfield = NA,
xaxis.prefix = "",
xaxis.unit = 1,
xaxis.suffix = "",
xaxis.round = 3,
xaxis.cex.label = 1,
xaxis.newline = TRUE,
xaxis.chronly = FALSE,
xaxis.width= TRUE,
xaxis.interval = 'auto',
xaxis.label.angle = 0,
xaxis.ticklen = 1,
xaxis.cex.tick = 1,
sep.lty = 2,
sep.lwd = 1,
sep.bg.col = 'gray95',
sep.draw = TRUE,
y0 = NA,
y1 = NA,
m.sep.lwd = 1,
# m.bg.col = NA, ## DEPRECATED
cmap.min = NA,
cmap.max = NA,
labels.suppress = FALSE,
labels.suppress.grl = labels.suppress,
labels.suppress.gr = labels.suppress,
bg.col = 'white', ## background color of whole thing
formatting = NA) {
## TODO: FIX THIS USING formals() and some eval / do.call syntax or something similar
new('gTrack', data = data, y.field = y.field, mdata = mdata, name = name, format = formatting,
edges = edges, vars = vars, draw.paths = draw.paths, colormaps = colormaps, height = height, ygap = ygap,
path.col = path.col,
path.lwd = path.lwd,
stack.gap = stack.gap, col = col, border = border, angle = angle, draw.var = draw.var,
gr.colorfield = gr.colorfield, y.quantile = y.quantile,
cex.label = cex.label, gr.cex.label.gr = gr.cex.label, gr.srt.label = gr.srt.label,
y.cap = y.cap, lwd.border = lwd.border, hadj.label = hadj.label, vadj.label = vadj.label, smooth = smooth,
round = round, ywid = ywid, ypad = ypad, seqinfo = seqinfo, circles = circles, lines = lines,
bars = bars, triangle = triangle, ylab = ylab, max.ranges = max.ranges, source.file.chrsub = source.file.chrsub,
y0.bar = y0.bar, yaxis = yaxis, yaxis.pretty = yaxis.pretty, yaxis.cex = yaxis.cex,
chr.sub = chr.sub, edgevars = edgevars, gr.labelfield = gr.labelfield,
grl.labelfield = grl.labelfield, grl.colorfield = gr.colorfield, grl.cexfield, xaxis.prefix = xaxis.prefix, xaxis.unit = xaxis.unit,
xaxis.suffix = xaxis.suffix, xaxis.round = xaxis.round, xaxis.interval = xaxis.interval,
xaxis.cex.label = xaxis.cex.label, xaxis.newline = xaxis.newline,
xaxis.chronly = xaxis.chronly, xaxis.width = xaxis.width,
labels.suppress = labels.suppress, labels.suppress.gr = labels.suppress.gr, labels.suppress.grl = labels.suppress.grl,
xaxis.label.angle = xaxis.label.angle, xaxis.ticklen = xaxis.ticklen,
xaxis.cex.tick = xaxis.cex.tick, sep.lty = sep.lty, sep.lwd = sep.lwd, sep.bg.col = sep.bg.col,
sep.draw = sep.draw, y0 = y0, y1 = y1, m.sep.lwd = m.sep.lwd, #m.bg.col = m.bg.col,
cmap.min = cmap.min, cmap.max = cmap.max, bg.col = bg.col)
}
setValidity('gTrack', function(object)
{
problems = c();
ALLOWED.DATA.CLASSES = c('GRangesList', 'GRanges', 'CompressedGRangesList', 'RleList', 'SimpleRleList', 'character', 'ffTrack');
if (length(object@data)>0)
{
if (!all(sapply(object@data, function(x) class(x) %in% ALLOWED.DATA.CLASSES)))
problems = c(problems, (paste('One or more of the provided data objects in data slot are not of the allowed data classes', paste(ALLOWED.DATA.CLASSES, collapse = ", "))))
else
{
# for (i in 1:nrow(object@formatting))
# if (!(object@formatting$format[i] %in% ALLOWED.FORMATS[[class(object@data[[i]])]]))
# problems = c(problems, sprintf('Format specified for data track %s (%s) is not allowed. Available formats for this track are: %s.', i, object@formatting$format[i], paste(ALLOWED.FORMATS[[class(object@data[[i]])]], collapse = ", ")))
if (length(object@data) != nrow(object@formatting))
problems = c(problems, 'Length of data is not equal to the number of formatting entries')
}
}
else
{
if (nrow(object@formatting) != 0)
problems = c(problems, 'Null trackdata object has incompatible fields');
}
if (!all(sapply(object@edges, is.data.frame)))
problems = c(problems, 'Some trackdata edges attributes are not data.frames')
else if (any(!sapply(object@edges, function(x) if (nrow(x)>0) all(c('from', 'to') %in% colnames(x)) else T)))
problems = c(problems, 'Some nonempty trackdata edges attributes are missing $to and $from fields')
else if (any(!sapply(1:length(object@data), function(x)
if (nrow(object@edges[[x]])>0)
all(object@edges[[x]]$from <= length(object@data[[x]])) & all(object@edges[[x]]$to <= length(object@data[[x]]))
else T)))
problems = c(problems, 'Some nonempty trackdata edges $to and $from fields are out of bounds (ie exceed the length of the data field of the corresponding gTrack item')
if (!is.null(formatting(object)$y.field) && !all(is.na(formatting(object)$y.field)))
{
nix = !is.na(object$y.field) & sapply(dat(object), inherits, 'GRanges')
if (any(nix))
tmp = sapply(which(nix), function(x) object$y.field[x] %in% names(values(dat(object)[[x]])))
else
tmp = TRUE
if (any(!tmp))
{
problems = c(problems, paste('These y.fields are not found in their respective GRanges object:',
paste(object$y.field[nix][!tmp], collapse = ', ')))
}
}
if (length(object@data) != length(object@colormap))
problems = c(problems, 'Length of object is not the same length as colormap')
if (any(ix <- sapply(object@data, class) == 'character'))
{
for (iix in which(ix))
{
x = object@data[[iix]]
if (grepl('(\\.bw$)|(\\.bigwig$)', x, ignore.case = T))
{
f = BigWigFile(normalizePath(x))
slen = tryCatch(GenomeInfoDb::seqlengths(f), error = function(x) NULL)
if (is.null(slen))
problems = c(problems, sprintf('External file %s is not a valid and existing .bw / .bigwig file', x))
}
else if (grepl('\\.wig', x, ignore.case = T))
{
f = WIGFile(normalizePath(x))
slen = tryCatch(GenomeInfoDb::seqlengths(f), error = function(x) NULL)
if (is.null(slen))
problems = c(problems, sprintf('External file %s is not a valid and existing .wig file', x))
}
else if (grepl('\\.bed', x, ignore.case = T))
{
f = BEDFile(normalizePath(x))
# slen = tryCatch(GenomeInfoDb::seqlengths(f), error = function(x) NULL)
# if (is.null(slen))
# problems = c(problems, sprintf('External file %s is not a valid and existing .bed file', x))
}
else if (grepl('\\.gff', x, ignore.case = T))
{
f = GFFFile(normalizePath(x))
slen = tryCatch(GenomeInfoDb::seqlengths(f), error = function(x) NULL)
if (is.null(slen))
problems = c(problems, sprintf('External file %s is not a valid and existing .gff file', x))
}
else if (grepl('\\.2bit', x, ignore.case = T))
{
f = TwoBitFile(normalizePath(x))
slen = tryCatch(GenomeInfoDb::seqlengths(f), error = function(x) NULL)
if (is.null(slen))
problems = c(problems, sprintf('External file %s is not a valid and existing .2bit file', x))
}
else if (grepl('\\.bedgraph', x, ignore.case = T))
{
f = BEDGraphFile(normalizePath(x))
# slen = tryCatch(GenomeInfoDb::seqlengths(f), error = function(x) NULL)
# if (is.null(slen))
# problems = c(problems, sprintf('External file %s is not a valid and existing .bedgraph file', x))
}
else
problems = c(problems, sprintf('External file %s does not map to a supported UCSC format or .rds. Supported files must have one of the following extensions: .bw, .bed. .bedgraph, .2bit, .wig, .gff, .rds.', x))
}
}
if (length(problems)==0)
TRUE
else
problems
}
);
#suppressWarnings(removeMethod('[', 'gTrack')) ## takes care of stupid R 2.15 bug
#' @name [
#' @title [
#' @description
#'
#' Subsetting tracks of gTrack gt
#' gt[1]
#' gt[ix] # where ix is an integer vector
#'
#' @param x \code{gTrack} object
#' @param i Integer specifying which \code{gTrack} to grab
#' @docType methods
#' @export
#' @aliases [,gTrack,ANY,ANY,ANY-method
#' @rdname sub-methods
#' @author Marcin Imielinski
setMethod('[', 'gTrack', function(x, i)
{
if (is.logical(i))
i = which(i)
x@data = x@data[i]
x@seqinfo = x@seqinfo
# x@height = x@height[i]
x@formatting = x@formatting[i, ]
x@colormap = x@colormap[i]
x@edges = x@edges[i]
if (.hasSlot(x, 'mdata'))
x@mdata = x@mdata[i]
else
x@mdata = rep(list(matrix()), length(x))
return(x)
})
#' @name mdata
#' @title mdata
#' @description
#'
#' Accessing submatrices of mdata object for length 1 gTrack gt
#'
#' Usage: (igr, jgr are GRanges objects corresponding to slices of matrix to be accessed from gt)
#' mdata(gt, igr, jgr)
#' @param x \code{gTrack} object containing a matrix in the \code{mdata} slot
#' @param igr \code{GRegion} object specifying positions to subset matrix with
#' @param jgr \code{GRegion} object specifying positions to subset matrix with
#' @author Jeremiah Wala
#' @docType methods
#' @rdname mdata-methods
#' @export
setGeneric('mdata', function(x, igr=NULL, jgr=NULL) standardGeneric('mdata'))
#' @rdname mdata-methods
#' @aliases mdata,gTrack,ANY,ANY,ANY-method
setMethod('mdata', signature=c("gTrack", "ANY", "ANY"), function(x, igr = NULL, jgr = igr)
{
if (is.null(igr) & is.null(jgr))
return(x@mdata)
if (is.null(igr))
igr = si2gr(x)
if (is.null(jgr))
jgr = si2gr(x)
if (is(igr, 'Rle'))
igr = as.character(igr)
if (is(jgr, 'Rle'))
jgr = as.character(jgr)
if (is.character(igr))
igr = si2gr(x)[igr]
if (is.character(jgr))
jgr = si2gr(x)[jgr]
if (length(igr)==0 | length(jgr)==0)
return(NULL)
out = lapply(1:length(x), function(y)
{
if (is.null(x@mdata[[y]]))
return(NULL)
i = gUtils::gr.in(x@data[[y]], igr)
j = gUtils::gr.in(x@data[[y]], jgr)
rown = dedup(gr.string(x@data[[y]][i]))
coln = dedup(gr.string(x@data[[y]][j]))
tmp = (x@mdata[[y]][i, j] + t(x@mdata[[y]][j, i]))/2
rownames(tmp) = rown
colnames(tmp) = coln
return(tmp)
})
if (length(x)==1)
out = out[[1]]
return(out)
})
#' @name $
#' @title $
#' @description
#'
#' Accessing columns of gTrack formatting data.frame
#'
#' @param x \code{gTrack} object
#' @param name Name of the \code{formatting} field to view
#' @docType methods
#' @rdname cash-methods
#' @aliases $,gTrack-method
#' @export
#' @author Marcin Imielinski
setMethod('$', 'gTrack', function(x, name)
{
# return(x@formatting[, name])
return(x@formatting[[name]])
})
#' @name $<-
#' @title $<-
#' @description
#'
#' Setting formats of gTrack object - ie modifying the formatting(gt) data frame after
#' an object has already been instantiated
#'
#' Usage:
#' gt$y.field = 'score'
#' gt$gr.colorfield[1] = 'readtype'
#'
#' @param x \code{gTrack} object to alter \code{formatting} field of
#' @param name \code{formatting} field to alter
#' @param value New value
#' @docType methods
#' @rdname cash-set-methods
#' @aliases $<-,gTrack-method
#' @export
#' @author Marcin Imielinski
setMethod('$<-', 'gTrack', function(x, name, value)
{
x@formatting[, name] = value
return(x)
})
#' @name length
#' @title length
#' @description
#'
#' Getting length of gTrack object gt
#' Usage:
#' length(gt)
#' @param x \code{gTrack} object
#' @docType methods
#' @rdname length-methods
#' @aliases length,gTrack-method
#' @export
#' @author Marcin Imielinski
setMethod('length', 'gTrack', function(x)
{
return(length(x@data))
})
#' @name reduce
#' @title reduce
#' @description
#'
#' Computing the GRanges footprint of the gTrack object on the genome
#' usage:
#' reduce(gt) # outputs a GRanges
#' @param x \code{gTrack} object to retrieve reduced \code{GRanges} from
#' @param ... additional arguments to GRanges reduce function
#' @return \code{GRanges} with the minimal footprint of the \code{gTrack} data
#' @importFrom GenomicRanges reduce
#' @importFrom gUtils gr.sub
#' @docType methods
#' @rdname reduce-methods
#' @aliases reduce,gTrack-method
#' @export
#' @author Marcin Imielinski
setMethod('reduce', 'gTrack', function(x, ... )
{
if (length(x)==0)
return(GRanges(seqlengths = GenomeInfoDb::seqlengths(x)))
.dat2gr = function(y)
{
if (is(y, 'Rle'))
{
y[is.na(y)] = 0;
out = as(y!=0, 'GRanges');
return(out[out$score])
}
else if (is(y, 'GRangesList'))
gr.stripstrand(unlist(y))
else
gr.stripstrand(y)
}
if (length(x)==1)
return(reduce(.dat2gr(x@data[[1]]), ...))
else
return(reduce(do.call('grbind', lapply(x@data, .dat2gr)), ... ))
})
#uppressWarnings(removeMethod('seqinfo', 'gTrack')) ## takes care of stupid R 2.15 bug
#' @name seqinfo
#' @title seqinfo
#' @description
#'
#' returns Seqinfo of gTrack object gt
#' Usage:
#' seqinfo(gt)
#' @docType methods
#' @param x \code{gTrack} object
#' @return \code{seqinfo}
#' @importFrom GenomicRanges seqinfo
#' @export
#' @author Marcin Imielinski
setMethod("seqinfo", signature(x = "gTrack"), function(x)
{
return(x@seqinfo)
})
#' @name seqinfo<-
#' @title seqinfo,-
#' @description
#'
#' set seqinfo property of gTrack
#'
#' @export
#' @param .Object \code{gTrack} object
#' @param value new \code{Seqinfo} object
#' @docType methods
#' @rdname seqinfo-set-methods
#' @author Marcin Imielinski
setGeneric('seqinfo<-', function(.Object, value) standardGeneric('seqinfo<-'))
#' @rdname seqinfo-set-methods
#' @aliases seqinfo,gTrack-method
setReplaceMethod('seqinfo', 'gTrack', function(.Object, value)
{
.Object@seqinfo = value;
.Object@data = lapply(dat(.Object), gr.fix, value)
validObject(.Object)
return(.Object)
});
#' @name c
#' @title c
#' @description
#'
#' Concatenate gTrack objects gt1, gt2, gt3
#' c(gt1, gt2, gt3) # returns a gTrack object with the component tracks "stacked"
#'
#' @param x Initial \code{gTrack} object
#' @param ... Any number of \code{gTrack} objects
#' @param recursive If recursive = TRUE, the function recursively descends through lists (and pairlists) combining all their elements into a vector [FALSE]
#' @docType methods
#' @rdname c-methods
#' @aliases c,gTrack-method
#' @export
#' @author Marcin Imielinski
setMethod('c', 'gTrack', function(x, ..., recursive = FALSE)
{
args = list(x, ...)
if (any(ix <- sapply(args, inherits, 'trackData')))
args[ix] = lapply(args[ix], function(y) {class(y) = 'gTrack'; return(y)})
if (any(!sapply(args, inherits, 'gTrack')))
stop('some objects in gTrack concatenation are not gTrack')
has.mdata = sapply(args, .hasSlot, 'mdata')
if (any(!has.mdata))
for (j in which(!has.mdata))
args[[j]]@mdata = rep(list(matrix()), length(args[[j]]))
args = args[!sapply(args, is.null)]
out <- gTrack(data = do.call('c', lapply(1:length(args), function(y) args[[y]]@data)),
# seqinfo = do.call('c', lapply(1:length(args), function(y) args[[y]]@seqinfo)),
colormaps = do.call('c', lapply(args, function(y) y@colormap)),
edges = do.call('c', lapply(args, function(y) y@edges)),
mdata = do.call('c', lapply(1:length(args), function(y) args[[y]]@mdata)),
format = do.call('rrbind', lapply(args, formatting)),
y.field = do.call('c', lapply(args, function(y) formatting(y)$y.field)),
yaxis = do.call('c', lapply(args, function(y) formatting(y)$yaxis)))
#out@mdata = do.call('c', lapply(1:length(args), function(y) args[[y]]@mdata))
if (is.null(out$name))
out$name = NA
## name vs track.name --> get rid of track.name
if (!is.null(out$track.name))
out$name = ifelse(is.na(out$name), out$track.name, out$name)
out$track.name = NULL
## formatting(out) <- do.call('rrbind', lapply(args, formatting))
if ('is.null' %in% names(formatting(out)))
{
is.n = out$is.null
is.n[is.na(is.n)] = F
out = out[!is.n]
}
return(out)
})
#' @name formatting
#' @title formatting
#' @description
#'
#' Get data frame specifying formatting of gTrack object gt
#' usage:
#' formatting(gt)
#' gt # will just display the formatting data.frame
#'
#' If you want to access particular fields of the formatting data.frame, just use the "$" accessor like you would for a data.frame
#'
#' @param .Object \code{gTrack} object to extracting the formatting data.frame from
#' @docType methods
#' @rdname formatting-methods
#' @export
#' @author Marcin Imielinski
setGeneric('formatting', function(.Object) standardGeneric('formatting'))
#' @name xaxis
#' @title Retrieves the xaxis parameters
#' @description
#'
#' Return the portion of the gTrack @format field responsible for
#' formatting the x-axis
#' @param .Object \code{gTrack} object to retrieve xaxis parameters from
#' @return \code{data.frame} which is a subset of \code{formatting} showing xaxis params
#' @author Jeremiah Wala
#' @docType methods
#' @rdname xaxis-methods
#' @export
setGeneric('xaxis', function(.Object) standardGeneric('xaxis'))
#' @rdname xaxis-methods
#' @aliases xaxis,gTrack-method
setMethod('xaxis', 'gTrack', function(.Object) {
xaxis.fields <- c("xaxis.prefix", "xaxis.suffix", "xaxis.unit",
"xaxis.round", "xaxis.interval", "xaxis.pos",
"xaxis.nticks", "xaxis.pos.label",
"xaxis.cex.label","xaxis.newline", "xaxis.chronly",
"xaxis.width", "xaxis.label.angle","xaxis.ticklen")
return(.Object@formatting[,which(colnames(.Object@formatting) %in% xaxis.fields)])
})
#' @name sep
#' @title Retrieves the seperator graphical parameters
#' @description
#'
#' Return the portion of the gTrack @format field responsible for
#' formatting the windows and their separators
#' @author Jeremiah Wala
#' @param .Object \code{gTrack} object to extract separator params from
#' @return \code{data.frame} with the seperator paramters
#' @docType methods
#' @rdname sep-methods
#' @export
setGeneric('sep', function(.Object) standardGeneric('sep'))
#' @rdname sep-methods
#' @aliases sep,gTrack-method
setMethod('sep', 'gTrack', function(.Object) {
sep.fields <- c("sep.lty", "sep.lwd", "sep.bg.col", "sep.draw")
return(.Object@formatting[,which(colnames(.Object@formatting) %in% sep.fields)])
})
#' @rdname formatting-methods
#' @aliases formatting,gTrack-method
setMethod('formatting', 'gTrack', function(.Object)
{
return(.Object@formatting)
})
#' @name edgs
#' @title edgs
#' @description
#'
#' Get edges list associated with gTrack object, is is a list of data.frames
#' each which contains a data.frame specifying formatted edges, as pairs of indices into the
#' underlying GRanges object $from, $to and additional formats $col $lwd $lty
#'
#' usage:
#' edgs(gt)
#' @param .Object \code{gTrack} object to extract the edge list from
#' @docType methods
#' @rdname edgs-methods
#' @export
setGeneric('edgs', function(.Object) standardGeneric('edgs'))
#' @rdname edgs-methods
#' @aliases edgs,gTrack-method
setMethod('edgs', 'gTrack', function(.Object)
{
# .Object = list(...)[[1]]
return(.Object@edges)
})
setGeneric('vars', function(.Object) standardGeneric('vars'))
#' @name vars
#' @title vars
#' @description
#'
#' Get variants associated with track
#'
#' @keywords internal
#' @author Marcin Imielinski
setMethod('vars', 'gTrack', function(.Object)
{
# .Object = list(...)[[1]]
return(.Object@vars)
})
#' @name edgs<-
#' @title edgs<-
#' @description
#' Set edges data.frame associated with associated with gTrack object, the data.frame
#' that is being used to replace must have fields $from, $to, and can have optional fields $lwd
#' $lty, $col specifyign color and line type.
#'
#' usage:
#' edgs(gt)[[1]] <- new.edges.
#'
#' @param .Object \code{gTrack} object on which to set new edges
#' @param value New value of the edges
#' @docType methods
#' @rdname edgs-set-methods
#' @export
#' @author Marcin Imielinski
setGeneric('edgs<-', function(.Object, value) standardGeneric('edgs<-'))
#' @rdname edgs-set-methods
#' @aliases edgs<-,gTrack-method
setMethod('edgs<-', 'gTrack', function(.Object, value)
{
if (!all(sapply(value, is.data.frame)))
stop('edges attribute must be a list of data.frames')
non.empty = sapply(value, nrow)>0
if (any(non.empty) )
if (!all(sapply(value[non.empty], function(x) all(c('from', 'to') %in% colnames(x)))))
stop('data.frames of edges field of gTrack must be either empty or have $from and $to fields specified')
.Object@edges = value
validObject(.Object)
return(.Object)
})
#' @name clear
#' @title clear
#' @description
#'
#' Clear data from gTrack object (makes into an empty container with the same seqinfo and display features)
#' usage:
#' clear(gTrack)
#'
#' @param .Object \code{gTrack} object to clear
#' @author Marcin Imielinski
#' @docType methods
#' @rdname clear-methods
#' @export
setGeneric('clear', function(.Object) standardGeneric('clear'))
#' @rdname clear-methods
#' @aliases clear,gTrack-method
setMethod('clear', 'gTrack', function(.Object)
{
.Object = .Object[1]
.Object@data[[1]] = .Object@data[[1]][NULL]
validObject(.Object)
return(.Object);
})
#' @name dat
#' @title dat
#' @description
#'
#' Extract list of data contained inside a gTrack object gt. Each list item will contain
#' either a GRanges, RleList, or path to a file.
#'
#' usage
#' dat(gt)
#'
#' @param .Object \code{gTrack} to retrive data from
#' @docType methods
#' @rdname dat-methods
#' @export
#' @author Marcin Imielinski
setGeneric('dat', function(.Object) standardGeneric('dat'))
#' @rdname dat-methods
#' @aliases dat,gTrack-method
setMethod('dat', 'gTrack', function(.Object)
{
return(.Object@data)
})
#' @name formatting<-
#' @title formatting<-
#' @description
#'
#' Set formatting of gTrack object gt
#'
#' usage:
#' #gt is a gTrack object
#' formatting(td)$height <- 2
#'
#' #can also just use $ directly, like for a data frame, which is more convenient
#' td$height <- 2
#' @param .Object \code{gTrack} to set the formatting field for
#' @param value \code{data.frame} with the new formatting information
#' @docType methods
#' @rdname formatting-set-methods
#' @export
setGeneric('formatting<-', function(.Object, value) standardGeneric('formatting<-'))
#' @rdname formatting-set-methods
#' @aliases formatting<-,gTrack-method
setReplaceMethod('formatting', 'gTrack', function(.Object, value)
{
REQUIRED.COLUMNS = c('height', 'col', 'ygap', 'y.field');
if (nrow(value) != length(.Object))
stop('Replacement data frame has %s rows and the object has %s', nrow(value), length(value))
if (length(setdiff(REQUIRED.COLUMNS, colnames(value))))
stop('Replacement data frame is missing some columns')
.Object@formatting = value;
validObject(.Object)
return(.Object)
});
#' @name colormap
#' @title colormap
#' @description
#' Access colormap of gTrack object, this is a named list of named character vectors that specifies
#' the field of the underlying GRanges object that will be used to map a set of values
#' to a set of colors.
#' usage:
#' colormap(gt)
#'
#' @param .Object \code{gTrack} object to retrieve colormap from
#' @docType methods
#' @rdname colormap-methods
#' @export
#' @author Marcin Imielinski
setGeneric('colormap', function(.Object) standardGeneric('colormap'))
#' @rdname colormap-methods
#' @aliases colormap,gTrack-method
setMethod('colormap', 'gTrack', function(.Object)
{
return(.Object@colormap)
})
#' @name colormap<-
#' @title colormap<-
#' @description
#' Set colormap of gTrack object, this is a named list of named character vectors that specifies
#' the field of the underlying GRanges object that will be used to map a set of values
#' to a set of colors.
#'
#' usage:
#' colormap(gt)[1] = list(tumortype = c(lung = 'red', pancreatic = 'blue', colon = 'purple'))
#' @param .Object \code{gTrack} to set the colormap for
#' @param value New \code{colormap}
#' @docType methods
#' @rdname colormap-set-methods
#' @export
setGeneric('colormap<-', function(.Object, value) standardGeneric('colormap<-'))
#' @rdname colormap-set-methods
#' @aliases colormap<-,gTrack-method
setReplaceMethod('colormap', 'gTrack', function(.Object, value)
{
.Object@colormap = value;
validObject(.Object)
return(.Object)
});
#' @name show
#' @title show
#' @description Display a \code{gTrack} object
#' @docType methods
#' @param object \code{gTrack} to display
#' @rdname show-methods
#' @aliases show,gTrack-method
#' @export
#' @author Marcin Imielinski
setMethod('show', 'gTrack', function(object)
{
cat(sprintf('gTrack object with %s tracks with formatting:\n', length(object)))
print(formatting(object))
})
### utility function to allow "softer" matching of seqinfo's of genomes in chains
.identical.seqinfo = function(a, b)
{
df.a = data.frame(seqnames = GenomeInfoDb::seqnames(a), seqlengths = GenomeInfoDb::seqlengths(a), stringsAsFactors = F);
df.b = data.frame(seqnames = GenomeInfoDb::seqnames(b), seqlengths = GenomeInfoDb::seqlengths(b), stringsAsFactors = F);
df.a = df.a[order(df.a$seqnames), ]
df.b = df.b[order(df.b$seqnames), ]
rownames(df.a) = NULL
rownames(df.b) = NULL
return(identical(df.a, df.b))
}
## if (!isGeneric("plot"))
## setGeneric("plot", function(x, ...) standardGeneric("plot"))
#' Improved \code{rbind} for intersecting/union columns of \code{data.frames} or \code{data.tables}
#'
#' Like \code{rbind}, but takes the intersecting columns of the data.
#' @param ... Any number of \code{data.frame} or \code{data.table} objects
#' @param union Take union of columns (and put NA's for columns of df1 not in df2 and vice versa). \code{[TRUE]}
#' @param as.data.table Return the binded data as a \code{data.table}. \code{[FALSE]}
#' @return \code{data.frame} or \code{data.table} of the \code{rbind} operation
#' @importFrom data.table data.table rbindlist
#' @author Marcin Imielinski
rrbind = function (..., union = TRUE, as.data.table = FALSE)
{
dfs = list(...)
dfs = dfs[!sapply(dfs, is.null)]
dfs = dfs[sapply(dfs, ncol) > 0]
if (any(mix <- sapply(dfs, class) == 'matrix'))
dfs[mix] = lapply(dfs, as.data.frame)
names.list = lapply(dfs, names)
cols = unique(unlist(names.list))
unshared = lapply(names.list, function(x) setdiff(cols, x))
ix = which(sapply(dfs, nrow) > 0)
if (any(sapply(unshared, length) != 0))
expanded.dts <- lapply(ix, function(x) {
tmp = dfs[[x]]
if (is.data.table(dfs[[x]]))
tmp = as.data.frame(tmp)
tmp[, unshared[[x]]] = NA
return(data.table::as.data.table(as.data.frame(tmp[,
cols])))
})
else expanded.dts <- lapply(dfs, function(x) as.data.table(as.data.frame(x)[, cols]))
rout <- tryCatch(rbindlist(expanded.dts), error = function(e) NULL)
if (is.null(rout))
rout = data.table::as.data.table(do.call("rbind", lapply(expanded.dts,
as.data.frame)))
if (!as.data.table)
rout = as.data.frame(rout)
if (!union) {
shared = setdiff(cols, unique(unlist(unshared)))
rout = rout[, shared]
}
return(rout)
}
## #' @name plot
## #' @title plot method
## #'
## #' @export
## plot = function(x, y, ...) {
## UseMethod("plot")
## }
#' @name plot.gTrack
#' @title plot
#'
#' @description
#' Plot gTrack object in multi-track genome browser view. gTracks will appear as stacked interval
#' plots, line / bar / scatter plots, node-edge graphs, or triangle plots depending on the formatting
#' settings and additional data (eg mdata) provided. gTracks can be drawn across several non-contiguous
#' "windows" of the genome
#'
#' Additional argument "links" takes a GRangesList of signed interval pairs (ie each item is length 2)
#' representing genomic junctions as input. The junctions will be pointing right for + intervals and
#' left for - intervals. The meta data of the links GRangesList ($lwd $col $lty) can be set
#' to modulate the color, witwidth, and line style of the junction links .
#'
#' usage:
#' display(gt, win) # where win is a GRanges object
#' display(gt) # this will show the entire span of seqinfo(gt)
#' display(gt, '3:1e6-2e6') ## can use UCSC style strings
#' display(gt, GRanges(20, IRanges(20e6, 21e6)))
#' display(gt, '2', links = ra) # here, the entire chromosome 2 is shown, with rearrangements on top
#'
#'
#' @param windows GRanges specifying windows to view (can also be GRangesList), default is whole genome
#' @param links GRangesList of signed locus pairs specifying links to draw above the plot,
#' optional GRangesList values meta data specify formatting of individual links:
#' $label (text label), $col (line color), $lwd (line weight), $lty (line style),
#' $arrow (arrow flag), $col.arrow (arrow color), $v (vertical spline bulge), $w (horizontal spline bulge)
#' @param gap scalar numeric specifying gap between windows (only relevant if windows has length>1). Units of the gap are in genome coordinates
#' @param max.ranges scalar numeric > 0 specifying max number of ranges to draw in a window (via sampling). If specified, overrides gTrack max.ranges formatting feature.
#' @param ..., additional last-minute formatting changes to the gtrack can be entered here (eg col = 'blue')
#'
#' @docType methods
#' @author Marcin Imielinski, Jeremiah Wala
#'
#' @export plot.gTrack
#' @export
"plot.gTrack" = function(x, ##pplot (for easy search)
y,
windows = si2gr(seqinfo(x)), ## windows to plot can be Granges or GRangesList
links = NULL, ## GRangesList of pairs of signed locations,
gap = NULL, ## spacing betwen windows (in bp)
y.heights = NULL, # should be scalar or length(windows) if windows is a GRangesList
y.gaps = NULL, # relative heights of gaps below and above stacks (xaxes will be drawn here)
cex.xlabel = 1,
cex.ylabel = 1,
max.ranges = NA, # parameter for max ranges to draw on canvas in each track (overrides formatting)
links.feat = NULL, # links features override for links (must be nrow 1 or length(links) data frame
verbose=FALSE,
legend.params = list(),
... ## additional args to draw.grl OR last minute formatting changes to gTrack object
)
{
if (!missing(y))
windows = y
.Object = x
if (!missing(y))
windows = y
if (is(windows, 'numeric') | is(windows, 'integer'))
windows = as.character(windows)
if (is(windows, 'character'))
windows = unlist(parse.grl(windows, seqlengths(seqinfo(.Object))))
if (is(windows, 'GRangesList'))
windows = unlist(windows)
## make sure we have min legend data
if (!"xpos" %in% names(legend.params))
legend.params$xpos = 0
if (!"ypos" %in% names(legend.params))
legend.params$ypos = 1
if (!"plot" %in% names(legend.params))
legend.params$plot = TRUE
win.gap = gap ## PATCH: recasting some variable names
new.plot = TRUE
window.segs = list();
dotdot.args = list(...);
## parse the wind<ows into GRanges
windows = format_windows(windows, .Object)
windows = windows[width(windows)>0]
## if totally empty, plot blank and leave
if(!length(windows)) {
plot.blank(bg.col = bg.col)
return()
}
## make sure gTrack has all fields that are expected later
.Object <- prep_defaults_for_plotting(.Object)
if (is.null(formatting(.Object)$legend))
formatting(.Object)$legend = TRUE
else (any(is.na(formatting(.Object)$legend)))
formatting(.Object)$legend[is.na(formatting(.Object)$legend)] = TRUE
if (is.null(formatting(.Object)$legend.title))
formatting(.Object)$legend.title = NA
## name vs track.name saga
if (is.null(formatting(.Object)$track.name))
{
if (!is.null(formatting(.Object)$name))
formatting(.Object)$track.name = formatting(.Object)$name
else if (!is.null(formatting(.Object)$y.field))
formatting(.Object)$track.name = formatting(.Object)$y.field
else
formatting(.Object)$track.name = NA
}
has.colormap = sapply(colormap(.Object), length)>0
has.colorfield = !is.na(formatting(.Object)$gr.colorfield)
formatting(.Object)$legend = ifelse(is.na(formatting(.Object)$legend), formatting(.Object)$legend == TRUE & (has.colormap | has.colorfield), formatting(.Object)$legend)
numlegends = sum(formatting(.Object)$legend)
which.legend = which(formatting(.Object)$legend)
.Object$legend.title = ifelse(!is.na(.Object$legend.title), .Object$legend.title,
ifelse(has.colormap, names(colormap(.Object))[1:length(.Object)],
ifelse(has.colorfield, .Object$gr.colorfield,
ifelse(!is.na(.Object$track.name), .Object$track.name, ''))))
## add last minute formatting changes to gTrack
if (length(dotdot.args)>0)
for (f in intersect(names(dotdot.args), names(formatting(.Object))))
formatting(.Object)[, f] = dotdot.args[[f]]
## set the window gap .. so that windows don't collide
if (is.null(win.gap))
win.gap = sum(as.numeric(width(windows)))/30
## get the height of the stacks
if (is.null(y.heights) | length(y.heights) != length(windows))
##y.heights = rep(1, length(windows)) ## old from when we had windows as GRangesList
y.heights <- 1
## set the gaps between the gTracks
if (is.null(y.gaps) )
y.gaps = y.heights*0.8
else if (length(y.gaps) != length(windows))
y.gaps = rep(y.gaps[1], length(windows))
## ensure that we don't plot too much
if (!is.na(max.ranges))
formatting(.Object)$max.ranges = pmin(max.ranges, formatting(.Object)$max.ranges, na.rm = TRUE)
oth.ix = 1:(length(windows)-1);
top.gaps = 0.5*y.gaps
bottom.gaps = 0.5*y.gaps
if (length(windows)==1)
oth.ix = c()
ylim.stacks = data.frame(start = c(bottom.gaps[1], bottom.gaps[1] + cumsum(y.heights[oth.ix] + top.gaps[oth.ix] + bottom.gaps[oth.ix+1])),
end = cumsum(y.heights + top.gaps + bottom.gaps) - top.gaps)
oth.ix = 1:(length(.Object)-1);
if (length(.Object)==1)
oth.ix = c()
tmp.top.gaps = 0.5 * formatting(.Object)$ygap
tmp.bottom.gaps = 0.5 * formatting(.Object)$ygap
tmp.ylim.subplot = data.frame(start = c(tmp.bottom.gaps[1], tmp.bottom.gaps[1] +
cumsum(formatting(.Object)$height[oth.ix] + tmp.top.gaps[oth.ix] + tmp.bottom.gaps[oth.ix+1])),
end = cumsum(formatting(.Object)$height + tmp.top.gaps + tmp.bottom.gaps) - tmp.top.gaps)
ylim = c(0, max(ylim.stacks$end)+top.gaps[length(top.gaps)])
ylim.parent=ylim
window.ylims = data.frame(start = rep(NA, length(windows)), end = NA);
new.axis = TRUE;
this.windows = windows
#end(this.windows) <- end(this.windows) + 1 ## +1 added
#this.windows = gUtils::streduce(windows[[i]]) ##gr.stripstrand(GenomicRanges::trim(windows[[i]]))
##if (!inherits(this.windows, 'GRanges'))
## this.windows = gUtils::si2gr(this.windows)
i=1
this.ylim.subplot = tmp.ylim.subplot;
this.ylim.subplot$start = affine.map(pmin(1, this.ylim.subplot$start), ylim = unlist(ylim.stacks[i, c('start', 'end')]), xlim = c(0, 1))
this.ylim.subplot$end = affine.map(pmin(1, this.ylim.subplot$end), ylim = unlist(ylim.stacks[i, c('start', 'end')]), xlim = c(0, 1))
this.tmp.bottom.gap = tmp.bottom.gaps[1]*(ylim.stacks$end[i]-ylim.stacks$start[i])
this.xaxis.pos = this.ylim.subplot$start[1]-bottom.gaps[i]*0-this.tmp.bottom.gap
this.xaxis.pos.label = this.ylim.subplot$start[1]-5*bottom.gaps[i]/6-this.tmp.bottom.gap
ylim.stacks[i, 'xaxis.pos'] = this.xaxis.pos
## loop through the gTracks
for (j in 1:length(.Object))
{
par(xpd = NA);
cmap = colormap(.Object)[[j]];
cfield = names(colormap(.Object))[j]
if (is.na(cfield))
cfield = formatting(.Object)$gr.colorfield[j]
if (length(cmap)==0)
cmap = NA
## get the data into GRanges or GRangesList format
pre.filtered = FALSE;
if (.Object@formatting$triangle[j])
pre.filtered = TRUE
tt <- extract_data_from_tmp_dat(.Object, j, this.windows)
.Object = tt$o
tmp.dat = tt$t
this.windows = tt$w
## subsample if we need to for enforcing max.ranges
if (!is.na(formatting(.Object)$max.ranges[j]) && formatting(.Object)$max.ranges[j] > 0) {
tt <- enforce_max_ranges(.Object, pre.filtered, j, tmp.dat, this.windows)
tmp.dat = tt$t
pre.filtered = tt$p
}
## flag to tell us whether data is pre-filtered to window (ie in fftrack or rlelist)
## adjust y0 .bar
if (is.null((formatting(.Object)$y0.bar[j])) || is.na((formatting(.Object)$y0.bar[j])))
formatting(.Object)$y0.bar[j] = NA
## smooth the y.field data
if (!is.na(formatting(.Object)$y.field[j]) && is(tmp.dat, 'GRanges') && !is.na(formatting(.Object)$smooth[j]))
tmp.dat <- smooth_yfield(.Object, j, tmp.dat)
## fix y limits and apply log transform if needed
if (!is.na(formatting(.Object)$y.field[j]) && (is.na(formatting(.Object)$y0[j]) || is.na(formatting(.Object)$y1[j])))
.Object <- format_yfield_limits(.Object, j, tmp.dat, pre.filtered, this.windows)
# if (formatting(.Object[j])$format != 'ranges')
# stop("violated assumption. need to fix")
all.args = as.list(formatting(.Object[j]))
all.args = c(dotdot.args, all.args[!names(all.args) %in% names(dotdot.args)])
# all.args = list(
# col = formatting(.Object)$col[j],
# ywid = formatting(.Object)$ywid[j],
# border = formatting(.Object)$border[j],
# lwd.border = formatting(.Object)$lwd.border[j],
# adj.label = c(formatting(.Object)$hadj.label[j],
# formatting(.Object)$vadj.label[j]),
# gr.adj.label = c(0.5,
# formatting(.Object)$vadj.label[j]),
# angle = formatting(.Object)$angle[j],
# y.pad = formatting(.Object)$ypad[j],
# circles = formatting(.Object)$circles[j],
# lines = formatting(.Object)$lines[j],
# bars = formatting(.Object)$bars[j],
# y.grid.cex = formatting(.Object)$yaxis.cex[j],
# edges = edgs(.Object)[[j]],
# y0.bar = formatting(.Object)$y0.bar[j],
# stack.gap = formatting(.Object)$stack.gap[j])
# na.fields = names(formatting(.Object))[sapply(1:ncol(formatting(.Object)), function(field) is.na(formatting(.Object)[j, field]))]
# other.fields = setdiff(names(formatting(.Object)), c('name', 'height', 'ygap', 'stack.gap', 'lift', 'split', 'angle', 'format', 'lwd.border', 'source.file', 'source.file.chrsub', 'ypad', 'ywid', 'border', 'col', 'hadj.label', 'vadj.label', 'y.field', 'round', 'cex.ylabel', 'y.quantile', 'max.ranges', 'yaxis', 'yaxis.cex', 'is.null', 'yaxis.pretty', names(all.args))) ## remove na fields and anything else that might mess up draw.grl
#
# other.formats = structure(names = other.fields,
# lapply(other.fields, function(x) formatting(.Object)[j, x]))
# all.args[names(other.formats)] = other.formats;
#
# all.args = all.args[setdiff(names(all.args), setdiff(na.fields, c('col')))]
this.y.field = formatting(.Object)$y.field[j]
this.y.grid = NA;
if (is.na(this.y.field) | !(this.y.field %in% names(values(tmp.dat))))
this.y = this.ylim.subplot[j, ]
else
{
if (is.null(formatting(.Object)$log))
formatting(.Object)$log = NA
if (!is.na(formatting(.Object)$log[j]))
if (formatting(.Object)$log[j])
{
if (!is.null(tmp.dat$ywid))
tmp.dat$ywid = log10(tmp.dat$ywid)
values(tmp.dat)[, this.y.field] = log10(values(tmp.dat)[, this.y.field])
formatting(.Object)[j, 'y0'] = log10(formatting(.Object)[j, 'y0'])
formatting(.Object)[j, 'y1'] = log10(formatting(.Object)[j, 'y1'])
}
range.y = NULL;
if (all(c('y0', 'y1') %in% names(formatting(.Object))))
{
if (!is.na(formatting(.Object)[j, 'y0']) & !is.na(formatting(.Object)[j, 'y1']))
range.y = c(formatting(.Object)[j, 'y0'], formatting(.Object)[j, 'y1'])
else if (!is.na(formatting(.Object)[j, 'y0']) & is.na(formatting(.Object)[j, 'y1']))
range.y = c(formatting(.Object)[j, 'y0'], max(setdiff(values(tmp.dat)[, this.y.field], c(Inf, -Inf)), na.rm = T))
else if (is.na(formatting(.Object)[j, 'y0']) & !is.na(formatting(.Object)[j, 'y1']))
range.y = c(min(setdiff(values(tmp.dat)[, this.y.field], c(Inf, -Inf)), na.rm = T), formatting(.Object)[j, 'y1'])
}
if (length(tmp.dat)>0)
{
if (!is.null(tmp.dat$ywid)) ## remove any weird infinite ywids
{
if (any(ix <- is.infinite(values(tmp.dat)$ywid)))
values(tmp.dat)$ywid[ix] = NA
}
else
values(tmp.dat)$ywid = NA
if (is.null(range.y)) ## if y range is empty then pull from data
range.y = range(setdiff(values(tmp.dat)[, this.y.field], c(Inf, -Inf)), na.rm = T);
# range.y = range(setdiff(values(dat(.Object)[[j]])[, this.y.field], c(Inf, -Inf)), na.rm = T);
}
## however if there is a single data value, then we need to scale appropriately
if (diff(range.y)==0)
{
if (any(ix <- !is.na(values(tmp.dat)$ywid))) ## use ywid
range.y = range.y + 5*max(values(tmp.dat)$ywid[ix])*c(-1, 1)
else # otherwise use some arbitrary proportion around the value
range.y = range.y + abs(range.y)*0.2*c(-1, 1)
}
this.y.ticks = pretty(range.y, formatting(.Object)$yaxis.pretty[j])
if (is.null(formatting(.Object)$y.cap))
formatting(.Object)$y.cap = NA
if (!is.na(formatting(.Object)$y.cap[j])) ## cap values from top and bottom
this.y = affine.map(values(tmp.dat)[, this.y.field], ylim = unlist(this.ylim.subplot[j, ]), xlim = range(this.y.ticks), cap = formatting(.Object)$y.cap[j])
else
this.y = affine.map(values(tmp.dat)[, this.y.field], ylim = unlist(this.ylim.subplot[j, ]), xlim = range(this.y.ticks), cap = TRUE)
## scale ywid to axes
if (!is.na(formatting(.Object)$ywid[j]))
{
if (is.na(formatting(.Object)$y.field[j])) ## if no y.field then ignore ywid
{
all.args$ywid = formatting(.Object)$ywid[j] = NA
}
else ## scale it
{
all.args$ywid = formatting(.Object)$ywid[j] = formatting(.Object)$ywid[j] * abs(diff(unlist(this.ylim.subplot[j, ]))) / diff(range(this.y.ticks))
}
}
# this.y = affine.map(values(dat(.Object)[[j]])[, this.y.field], ylim = unlist(this.ylim.subplot[j, ]), xlim = range(this.y.ticks))
## if need, bump the range to include ybar base
# if (formatting(.Object)$y0.bar[j] < min(unlist(this.ylim.subplot[j, ])))
# this.ylim.subplot[j,'start'] <- formatting(.Object)$y0.bar[j]
if (is.na(formatting(.Object)$y0.bar[j]))
formatting(.Object)$y0.bar[j] = 0
all.args$y0.bar = affine.map(formatting(.Object)$y0.bar[j], ylim = unlist(this.ylim.subplot[j, ]), xlim = range(this.y.ticks))
if (!is.na(formatting(.Object)$y.field[j]))
{
# make pretty grid in range.y
this.y.grid = structure(affine.map(this.y.ticks, ylim = unlist(this.ylim.subplot[j, ]), xlim = range(this.y.ticks)), names = this.y.ticks)
if (!is.na(formatting(.Object)$log[j]))
if (formatting(.Object)$log[j])
names(this.y.grid) = signif(10^this.y.ticks)
}
}
if (is.null(.Object[j]$bars))
formatting(.Object)$bars[j] = FALSE
else if (is.na(.Object[j]$bars))
formatting(.Object)$bars[j] = FALSE
if (.Object[j]$bars && is.na(all.args$y0.bar))
all.args$y0.bar = this.ylim.subplot[j, 1]
if (is.null(.Object$chr.sub))
.Object$chr.sub = FALSE
if (is.na(.Object$chr.sub[j]))
.Object$chr.sub[j] = FALSE
if (.Object[j]$chr.sub)
tmp.windows = gr.sub(windows, 'chr', '')
else
tmp.windows = this.windows
## fix legend params
this.legend.params = legend.params
if (!formatting(.Object)$legend[j])
this.legend.params$plot = FALSE
else
{
this.legend.params$xpos = NA
if (!is.null(formatting(.Object)$legend.xpos))
this.legend.params$xpos = is.formatting(.Object)$legend.xpos[j]
if (!is.null(formatting(.Object)$legend.ypos))
this.legend.params$ypos = formatting(.Object)$legend.ypos[j]
jj = match(j, which.legend)
if (is.na(this.legend.params$xpos))
this.legend.params$xpos = seq(0, 1, length.out = numlegends)[jj]
this.legend.params$xpos = seq(0, 1, length.out = numlegends)[jj]
this.legend.params$xjust = c(0, 0.5, 1)[as.integer(cut(this.legend.params$xpos, c(-0.01, 0.2, 0.8, 1)))]
this.legend.params$title = .Object$legend.title[j]
}
## remove field "y" from tmp.dat if it exists
if (is(tmp.dat, 'GRangesList'))
{
values(tmp.dat)$y = NULL
}
else
{
if (!is.null(tmp.dat$y))
tmp.dat$y = NULL
}
main.args <- list(grl=tmp.dat,y = this.y, ylim = ylim,
xaxis.pos = this.xaxis.pos,xaxis.pos.label = this.xaxis.pos.label,
win.gap = win.gap[i],windows = tmp.windows,
new.plot = new.plot, new.axis = new.axis,
ylim.subplot = {if (all(is.na(this.y.grid))) unlist(this.ylim.subplot[j, ]) else NULL},
gr.colorfield = cfield,gr.colormap = cmap,
legend = formatting(.Object)$legend[j],
y.grid = {if (formatting(.Object)$yaxis[j]) this.y.grid else NA},
verbose=verbose,
ylim.parent=ylim.parent,mdata=.Object@mdata[[j]],
leg.params = this.legend.params,
adj.label = c(formatting(.Object)$hadj.label[j],
formatting(.Object)$vadj.label[j]),
gr.adj.label = c(0.5,
formatting(.Object)$vadj.label[j]),
y.pad = formatting(.Object)$ypad[j],
y.grid.cex = formatting(.Object)$yaxis.cex[j],
edges = edgs(.Object)[[j]])
all.args <- c(main.args, all.args[!names(all.args) %in% names(main.args)])
## clear out na.args to make default setting simpler downstream
na.args = sapply(all.args, function(x) if(is.vector(x)) all(is.na(x)) else FALSE)
if (any(na.args))
all.args = all.args[!na.args]
# main.args = c(main.args, all.args[setdiff(names(all.args), names(main.args))]);
#
# other.args = dotdot.args
# other.args = other.args[setdiff(names(other.args), names(main.args))]
## make empty plot
if (new.plot)
{
blank.main.args <- all.args
###blank.main.args = main.args;
###blank.main.args[[1]] = blank.main.args[[1]][c()]
#blank.main.args$y = list(start = min(this.ylim.subplot$start), end = max(this.ylim.subplot$end))
# if (any(is.na(.Object@formatting$triangle)) & any(.Object@formatting$triangle))
# blank.main.args$y = list(start = min(this.ylim.subplot$start[is.na(.Object@formatting$triangle)]), end = max(this.ylim.subplot$end[is.na(.Object@formatting$triangle)])) ## JEREMIAH
# else
blank.main.args$grl <- GRanges()
blank.main.args$y = list(start=min(this.ylim.subplot$start), end = max(this.ylim.subplot$end))
blank.main.args$triangle=FALSE
blank.main.args$sep.draw=FALSE
do.call('draw.grl', blank.main.args)
##do.call('draw.grl', c(blank.main.args, other.args))
all.args$new.plot = FALSE
all.args$new.axis = FALSE
}
new.plot = FALSE
new.axis = FALSE
##arrg <- c(main.args, other.args)
##form <- as.list(formatting(.Object[j]))
##arrg <- c(arrg, form[!names(form) %in% names(arrg)]) ## add in remaining args
window.segs = list();
## window.segs[[i]] = do.call('draw.grl', c(main.args, other.args))
if (formatting(.Object[j])$triangle)
{
all.args$sigma= all.args$smooth
window.segs[[i]] <- do.call('draw.triangle', all.args[names(all.args) %in% c("grl","y","mdata","ylim.parent","windows","win.gap","sigma", "max.ranges",
"cmap.min","cmap.max", "m.sep.lwd","legend","leg.params",
"islog","gr.colormap")])
} else {
window.segs[[i]] <- do.call('draw.grl', all.args)
}
this.tname = formatting(.Object[j])$track.name
if (!is.na(this.tname))
{
ylab.las = if (length(dotdot.args[["ylab.las"]]) == 0) 0 else dotdot.args[["ylab.las"]]
ylab.adj = if (length(dotdot.args[["ylab.adj"]]) == 0) c(0.5, 1) else dotdot.args[["ylab.adj"]]
if (ylab.las %in% c(0,3)) {
ylab.angle = -90
} else if (ylab.las %in% c(1,2)) {
ylab.angle = 0
}
this.cex.ylabel = ifelse(!is.null(formatting(.Object[j])$cex.ylabel), formatting(.Object[j])$cex.ylabel, cex.ylabel)
text(par('usr')[2], mean(unlist(this.ylim.subplot[j, c('start', 'end')])),
this.tname, srt = ylab.angle, adj = ylab.adj, cex = this.cex.ylabel)
}
}
if (is.null(links))
links = GenomicRanges::GRangesList()
if (length(links)>0) # draw rearrangement links
{
if (!inherits(links, 'GRangesList'))
stop('links must be a GRangesList')
## make sure chr.sub works with objects
chr.sub = any(formatting(.Object)[,'chr.sub'], na.rm = TRUE)
if (!is.null(chr.sub) && !is.na(chr.sub) && chr.sub) {
links = gUtils::gr.sub(links, "chr", "")
}
# first map rearrangements to various windows>
win.u = this.windows
win.u$grl.ix = 1 ##holdover from grangeslist windows
##win.u = gr.stripstrand(grl.unlist(windows))
window.segs.u = do.call(rbind, window.segs)
window.segs.u$width = window.segs.u$end - window.segs.u$start + 1
window.segs.xlim = do.call('rbind', lapply(window.segs, function(x) data.frame(start = min(x$start), end = max(x$end))))
links.u = grl.unlist(links)
if (any(table(links.u$grl.ix)!=2))
stop('Links should be GRangesList of range pairs.')
links.p = grl.pivot(links)
## flip strands to conform to connectors convention (- connection to left side and + is connection to right side)
## with ra specification convention (- refers to segment to left of breakpoint and + refers to segment to right)
GenomicRanges::strand(links.p[[1]]) = c("*"= "-", "-" = "+", "+" = "-")[as.character(GenomicRanges::strand(links.p[[1]]))]
GenomicRanges::strand(links.p[[2]]) = c("*" = "+", "-" = "+", "+" = "-")[as.character(GenomicRanges::strand(links.p[[2]]))]
## find overlaps with windows and calculate their window specific coordinates
l1 = gr.findoverlaps(links.p[[1]], win.u)
values(l1) = cbind(as.data.frame(values(l1)), as.data.frame(values(links)[l1$query.id, , drop = FALSE]))
GenomicRanges::strand(l1) = GenomicRanges::strand(links.p[[1]])[l1$query.id]
l1$stack.id = win.u$grl.ix[l1$subject.id]
l1$y.pos = ylim.stacks$xaxis.pos[l1$stack.id]
l1$x.pos = mapply(function(x,y,z,a) (y-z)*a + x, x = window.segs.u[l1$subject.id,]$start, y = start(l1),
z = start(win.u[l1$subject.id]), a = window.segs.u$width[l1$subject.id] / width(win.u)[l1$subject.id])
l2 = gr.findoverlaps(links.p[[2]], win.u)
values(l2) = cbind(as.data.frame(values(l2)), as.data.frame(values(links)[l2$query.id, , drop = FALSE]))
GenomicRanges::strand(l2) = GenomicRanges::strand(links.p[[2]])[l2$query.id]
l2$stack.id = win.u$grl.ix[l2$subject.id]
l2$y.pos = ylim.stacks$xaxis.pos[l2$stack.id]
l2$x.pos = mapply(function(x,y,z,a) (y-z)*a + x, x = window.segs.u[l2$subject.id,]$start, y = start(l2),
z = start(win.u[l2$subject.id]), a = window.segs.u$width[l2$subject.id] / width(win.u)[l2$subject.id])
.fix.l = function(ll)
{
if (!is.null(links.feat))
for (col in names(links.feat))
values(ll)[, col] = links.feat[ll$query.id, col]
# set up connectors
if (is.null(ll$v))
ll$v = y.gaps[ll$stack.id]/4
else
ll$v = y.gaps[ll$stack.id]*ll$v/2
if (is.null(ll$h))
ll$h = (window.segs.xlim$end[ll$stack.id] - window.segs.xlim$start[ll$stack.id])/20
else
ll$h = (window.segs.xlim$end[ll$stack.id] - window.segs.xlim$start[ll$stack.id])*ll$h
if (is.null(ll$arrow))
ll$arrow = TRUE
if (is.null(ll$cex.arrow))
ll$cex.arrow = 1
if (is.null(ll$lwd))
ll$lwd = 1
if (is.null(ll$lty))
ll$lty = 3
if (is.null(ll$col))
ll$col = 'red'
if (is.null(ll$col.arrow))
ll$col.arrow = ll$col
return(ll)
}
if (length(l1)>0)
l1 = .fix.l(l1)
if (length(l2)>0)
l2 = .fix.l(l2)
## now pair up / merge l1 and l2 using query.id as primary key
if (length(l1)>0 & length(l2)>0)
{
pairs = merge(data.frame(l1.id = 1:length(l1), query.id = l1$query.id), data.frame(l2.id = 1:length(l2), query.id = l2$query.id))[, c('l1.id', 'l2.id')]
l1.paired = GenomicRanges::as.data.frame(l1)[pairs[,1], ]
l2.paired = GenomicRanges::as.data.frame(l2)[pairs[,2], ]
}
else
{
l2.paired = l1.paired = data.frame()
pairs = data.frame()
}
## some l1 and l2 will be unpaired
l.unpaired = GRanges(seqlengths = GenomeInfoDb::seqlengths(links));
p1 = p2 = c();
if (nrow(pairs)>0)
{
p1 = pairs[,1]
p2 = pairs[,2]
}
if (length(l1)>0)
l.unpaired = grbind(l.unpaired, l1[setdiff(1:length(l1), p1)])
if (length(l2)>0)
l.unpaired = grbind(l.unpaired, l2[setdiff(1:length(l2), p2)])
if (length(l.unpaired)>0)
{
l.unpaired$v = l.unpaired$v/4
l.unpaired$h = l.unpaired$h/2
l.unpaired$col = alpha(l.unpaired$col, 0.5)
## unpaired will get a "bridge to nowhere" in the proper direction (eg down)
l.unpaired$y.pos = ylim.stacks$end[l.unpaired$stack.id]
# l.unpaired$y.pos2 = l.unpaired$y.pos + top.gaps[l.unpaired$stack.id]
l.unpaired$y.pos2 = l.unpaired$y.pos + l.unpaired$v
connectors(l.unpaired$x.pos, l.unpaired$y.pos, as.character(GenomicRanges::strand(l.unpaired)),
l.unpaired$x.pos, l.unpaired$y.pos2,
as.character(GenomicRanges::strand(l.unpaired)),
v = abs(l.unpaired$v), h = l.unpaired$h, type = 'S',
f.arrow = l.unpaired$arrow, b.arrow = l.unpaired$arrow,
cex.arrow = 0.2*l.unpaired$cex.arrow,
col.arrow = l.unpaired$col.arrow,
lwd = l.unpaired$lwd, lty = l.unpaired$lty, col = l.unpaired$col)
if (!is.null(l.unpaired$label))
{
cex = 0.5;
if (!is.null(l.unpaired$cex.label))
cex = l.unpaired$cex.label
# text(l.unpaired$x.pos, l.unpaired$y.pos2, l.unpaired$label, adj = c(0.5, 0.5), cex = cex)
l.unpaired$text.y.pos =l.unpaired$y.pos - (ylim.stacks$end[l.unpaired$stack.id]-ylim.stacks$start[l.unpaired$stack.id])/100
text(l.unpaired$x.pos, l.unpaired$text.y.pos, l.unpaired$label, adj = c(0.5, 1), cex = cex)
}
}
## now draw connectors for paired links
## pairs on the same level will get a "U" link,
## pairs on different levels will get "S" links
## fix y distances so that connections go from topmost part of bottom most connection to the bottom most part of
## top connection (ie the xaxis position)
if (nrow(l1.paired)>0)
{
l1.paired$bottom = l1.paired$y.pos < l2.paired$y.pos
if (any(l1.paired$bottom))
l1.paired$y.pos[l1.paired$bottom] = ylim.stacks$end[l1.paired$stack.id[l1.paired$bottom]]
if (any(!l1.paired$bottom))
l2.paired$y.pos[!l1.paired$bottom] = ylim.stacks$end[l2.paired$stack.id[!l1.paired$bottom]]
# also make all c type connections top connections with positive v
ctype = ifelse(l1.paired$stack.id == l2.paired$stack.id, 'U', 'S')
l1.paired$y.pos[ctype == 'U'] = ylim.stacks$end[l1.paired$stack.id[ctype == 'U']]
l2.paired$y.pos[ctype == 'U'] = ylim.stacks$end[l2.paired$stack.id[ctype == 'U']]
l1.paired$v[ctype == 'U'] = abs(l1.paired$v[ctype == 'U'])
win.width = diff(par('usr')[1:2])
l1.paired$v = l1.paired$v * affine.map(abs(l2.paired$x.pos - l1.paired$x.pos)/diff(par('usr')[1:2]), xlim = c(0, 1), ylim = c(0.5, 1)) ## make longer links taller
connectors(l1.paired$x.pos, l1.paired$y.pos, l1.paired$strand, l2.paired$x.pos, l2.paired$y.pos, l2.paired$strand,
v = l1.paired$v, h = l1.paired$h, type = ctype,
f.arrow = l1.paired$arrow, b.arrow = l1.paired$arrow, cex.arrow = 0.2*l1.paired$cex.arrow, col.arrow = l1.paired$col.arrow,
lwd = l1.paired$lwd, lty = l1.paired$lty, col = l1.paired$col)
if (!is.null(l1.paired$label))
{
cex = 0.5;
if (!is.null(l1.paired$cex.label))
cex = l1.paired$cex.label
## text(l1.paired$x.pos, l1.paired$y.pos+l1.paired$v/2, l1.paired$label, adj = c(0.5, 0.5), cex = cex)
## text(l2.paired$x.pos, l2.paired$y.pos+l1.paired$v/2, l2.paired$label, adj = c(0.5, 0.5), cex = cex)
l1.paired$text.y.pos = l1.paired$y.pos - (ylim.stacks$end[l1.paired$stack.id]-ylim.stacks$start[l1.paired$stack.id])/100
l2.paired$text.y.pos = l2.paired$y.pos - (ylim.stacks$end[l2.paired$stack.id]-ylim.stacks$start[l2.paired$stack.id])/100
text(l1.paired$x.pos, l1.paired$text.y.pos, l1.paired$label, adj = c(0.5, 1), cex = cex)
text(l2.paired$x.pos, l2.paired$text.y.pos, l2.paired$label, adj = c(0.5, 1), cex = cex)
}
}
}
}
#})
###############
# Tracked data helper functions
#
###############
#' @name karyogram
#' @title karyogram
#' @description
#'
#' Returns gTrack displaying band pattern for hg18 or hg19 karyotype
#'
#' @param hg19 logical scalar flag, if T returns gTrack for hg19
#' @param bands logical scalar, if T returns gTrack with colored giemsa bands
#' @param arms logical scalar, if T and bands F, returns chromosome arms with different colors and centromeres and telomeres marked, if arms is F and bands F returns chromosomes, each with a different color
#' @param tel.width numeric scalar, specifies telomere width in bases (only relevant if arms = T, bands = F)
#' @param ... Additional arguments sent to the \code{gTrack} constructor
#' @export
#' @author Marcin Imielinski
karyogram = function(file = NULL, hg19 = TRUE, bands = TRUE, arms = TRUE, tel.width = 2e6, ... )
{
# ucsc.bands = get_ucsc_bands(hg19 = hg19)
if (is.null(file))
{
if (hg19)
ucsc.bands = readRDS(system.file("extdata", "ucsc.bands.hg19.rds", package = 'gTrack'))
else
ucsc.bands = readRDS(system.file("extdata", "ucsc.bands.hg18.rds", package = 'gTrack'))
}
else
{
ucsc.bands = dt2gr(fread(file, header = FALSE)[, .(seqnames = V1, start = V2, end = V3, band = V4, stain = V5)])
}
if (bands)
{
si = si2gr(ucsc.bands);
si = suppressWarnings(si[order(as.numeric(as.character(seqnames(si))))])
si = Seqinfo(as.character(seqnames(si)), width(si)+1)
ucsc.bands = sort(gr.fix(ucsc.bands, si, drop = T))
ucsc.bands = split(ucsc.bands, seqnames(ucsc.bands))
td = gTrack(list(ucsc.bands), colormaps = list(stain = c('gneg' = 'white', 'gpos25' = 'gray25', 'gpos50' = 'gray50', 'gpos75'= 'gray75', 'gpos100' = 'black', 'acen' = 'red', 'gvar' = 'pink', 'stalk' = 'blue')), border = 'black', ...)
formatting(td)$stack.gap[1] = 0
}
else
{
tmp = c(RColorBrewer::brewer.pal(11, 'BrBG'), RColorBrewer::brewer.pal(11, 'PiYG'), RColorBrewer::brewer.pal(11, 'RdYlGn'))
col.adjust = 10;
col.karyo = data.frame(
qtel = rep('black', length(tmp)),
qarm = lighten(tmp, col.adjust),
centro = rep('gray', length(tmp)),
parm = lighten(tmp, -col.adjust),
ptel = rep('black', length(tmp)), stringsAsFactors = F)
if (arms) ## draw arms with slightly different hues of the same color and black ranges for telomere / centromere
{
tmp.tel = aggregate(end ~ as.character(seqnames), data = GenomicRanges::as.data.frame(ucsc.bands), FUN = max)
tmp.tel = structure(tmp.tel[,2], names = as.character(tmp.tel[,1]))+1
telomeres = c(GRanges(names(tmp.tel), IRanges::IRanges(start = rep(1, length(tmp.tel)), end = rep(tel.width, length(tmp.tel))),
seqlengths = GenomeInfoDb::seqlengths(seqinfo(ucsc.bands))),
GRanges(names(tmp.tel), IRanges::IRanges(tmp.tel-tel.width+1, tmp.tel), seqlengths = GenomeInfoDb::seqlengths(seqinfo(ucsc.bands))))
values(telomeres)$lwd.border = 1;
values(telomeres)$border = 'black';
centromeres = reduce(ucsc.bands[values(ucsc.bands)$stain=='acen'])
values(centromeres)$lwd.border = 1;
values(centromeres)$border = 'black';
# arms = reduce(setdiff(ucsc.bands, centromeres))
arms = IRanges::gaps(c(telomeres, centromeres))
arms = arms[GenomicRanges::strand(arms)=='*']
values(arms)$lwd.border = 1;
values(arms)$border = 'black';
karyotype = sort(c(arms, centromeres, telomeres))
values(karyotype)$region = paste(as.character(seqnames(karyotype)), c('ptel', 'p', ' cen', 'q', 'qtel'), sep = '')
colmap = structure(as.vector(t(as.matrix(col.karyo)))[1:length(karyotype)], names = values(karyotype)$region)
GenomicRanges::strand(karyotype) = '+'
# karyotype = split(karyotype, seqnames(karyotype))
si = si2gr(karyotype);
si = suppressWarnings(si[order(as.numeric(as.character(seqnames(si))))])
si = Seqinfo(as.character(seqnames(si)), width(si)+1)
karyotype = gr.fix(karyotype, si, drop = T)
names(karyotype) = NULL;
td = gTrack(list(karyotype), colormaps = list(region= colmap), border = 'black',
stack.gap = 0, xaxis.interval = 1e7, ...)
}
else
{
karyotype = sort(reduce(ucsc.bands))
values(karyotype)$region = as.character(seqnames(karyotype))
colmap = structure(col.karyo$qarm[1:length(karyotype)], names = as.character(seqnames(karyotype)))
GenomicRanges::strand(karyotype) = '+'
si = si2gr(karyotype);
si = suppressWarnings(si[order(as.numeric(as.character(seqnames(si))))])
si = Seqinfo(as.character(seqnames(si)), width(si)+1)
karyotype = gr.fix(karyotype, si, drop = T)
names(karyotype) = NULL;
td = gTrack(list(karyotype), colormaps = list(region = colmap),
border = 'black', stack.gap = 0, xaxis.interval = 1e7, ...)
}
}
formatting(td)$angle = 10
# formatting(td)$ywid = 2
return(td)
}
#' @name file.url.exists
#' @title Check if a file or url exists
#' @param f File or url
#' @return TRUE or FALSE
#' @importFrom RCurl url.exists
#' @noRd
file.url.exists <- function(f) {
return(file.exists(f) || RCurl::url.exists(f))
}
#' @name read.rds.url
#' @title Checks if path is URL or file, then reads RDS
#' @param f File or url
#' @return data
#' @noRd
read.rds.url <- function(f) {
if (grepl("^http",f))
return(readRDS(gzcon(url(f))))
return(readRDS(f))
}
#' @name track.gencode
#' @title Constructor a \code{gTrack} from GENCODE transcripts
#'
#' @description
#' Returns gTrack object representing GENCODE transcripts and their components (utr, cds etc) with assigned colors.
#' Usually built from cached data objects but can also be built from provided GRangesList
#'
#' @param gencode (optional) path to GENCODE gtf or a GRanges object from which to generate and cache new gTrack GENCODE object (into path specified by GENCODE_DIR or provided as cached.dir argument to this function object). Only necessary to run if gTrack GENCODDE .rds has not already been cacned, (note: cached = FALSE in order for new GENCODE gTrack to be generated )
#' @param rg (optional) GRangesList representing transcript models imported from GENCODE gff3 file using rtracklayer import
#' @param genes (optional) character vector specifying genes to limit gTrack object to
#' @param gene.collapse scalar logical specifying whether to collapse genes by transcript (or used stored version of transcripts)
#' @param grep character vector for which to grep genes to positively select
#' @param build character can be 'hg19', 'hg38', if GENCODE_DIR and gencode not specified will download this build from mskilab.com
#' @param grepe character vector for which to grep genes which to exclude
#' @param bg.col scalar character representing background color for genespan
#' @param cds.col scalar character representing background color for CDS
#' @param cds.utr scalar character representing background color for UTR
#' @param st.col scalar character representing color of CDS start
#' @param en.col scalar character representing color of CDS end
#' @param cached logical scalar whether to use "cached" version provided with package
#' @param gr.srt.label scalar numeric specifying angle on exon label
#' @param gr.cex.label scalar numeric > 0 specifying character expansion on exon label
#' @param labels.suppress.gr scalar logical specifying whether to suppress exon label plotting
#' @param stack.gap stack.gap argument to gTrack
#' @param gencode.cols character vector specifying additional columns to include in the gencode data used to generate the gTrack. By default if gencode includes a metadata column "gene_label" then this column is included. "gene_label" is regarded as a special case in which the value for each gene is assumes to be unique and the "gene_label" value is also used to annotate the GRangesList associated with the output gTrack. This allows the users to supply alternative values for annotations of the genes. In this case gngt$grl.labelfield = 'gene_label' could be used in order for the "gene_label" values to be show up when plotting the gencode gTrack.
#' @param use.gene.ids (logical) use the gene_id column to name items in the GRangesList. By default set to FALSE to maintain backwards compatibility, which means "gene_name" is used to name items in the output. The problem with this is that gene_name is not unique. For example, 5S_rRNA is a gene_name that is associated with multiple different gene_id's, each representing a different copy of 5S_rRNA on various chromossomes.
#' @param ... additional arguments passed down to gTrack
#'
#' @export
#' @author Marcin Imielinski
track.gencode = function(gencode = NULL,
gene.collapse = TRUE,
genes = NULL,
grep = NULL,
grepe = NULL, ## which to exclude
build = 'hg19',
bg.col = alpha('blue', 0.1), cds.col = alpha('blue', 0.6), utr.col = alpha('purple', 0.4),
st.col = 'green',
en.col = 'red',
grl.labelfield, ## Don't touch these
gr.labelfield,
col,
cached = T, ## if true will use cached version
cached.dir = Sys.getenv('GENCODE_DIR'),
cached.path = paste(cached.dir, "gencode.composite.rds", sep = '/'), ## location of cached copy
cached.path.collapsed = paste(cached.dir, "gencode.composite.collapsed.rds", sep = '/'), ## location of cached copy
gr.srt.label = 0,
gr.cex.label = 0.3,
cex.label = 0.5,
labels.suppress.gr = T,
drop.rp11 = TRUE,
stack.gap = 1e6,
gencode.cols = 'gene_label',
use.gene.ids = FALSE,
...)
{
if (nchar(cached.dir)==0)
{
url <- sprintf("http://mskilab.com/gTrack/%s", build)
cached.path <- file.path(url, "gencode.composite.rds")
cached.path.collapsed <- file.path(url, "gencode.composite.collapsed.rds")
todl = ifelse(gene.collapse, cached.path.collapsed, cached.path)
if (!file.url.exists(todl))
stop(paste0("Can't successfully download file from ", todl, ". Please check if this URL exists or if there is an issue with your internet connection"))
warning('Downloading GENCODE track for genome build ', build, ' from ', ifelse(gene.collapse, cached.path.collapsed, url), ' - for quicker loading, download this file locally and point GENCODE_DIR environment variable to the enclosing directory. To generate and use a new cached gTrack object from a GENCODE .gtf or .gff3 set GENCODE_DIR env variable to an existing local directory and run track.gencode(url_or_path_to_gencode_gtf) once, and then use track.gencode() subsequently to access that cached gTrack object.')
}
warning(sprintf(paste(ifelse(file.url.exists(cached.path), "Pulling", "Caching"),
"gencode annotations from %s"),
ifelse(gene.collapse, cached.path.collapsed, cached.path)))
if (!cached | (!gene.collapse & !file.url.exists(cached.path)) | (gene.collapse & !file.url.exists(cached.path.collapsed))) ## if no composite refgene copy, then make from scratch
{
if (!file.exists(cached.dir))
{
message('Cached dir ', cached.dir, ' does not exist: creating ...')
system(paste('mkdir -p', cached.dir))
}
if (is.null(gencode))
gencode = readRDS(gzcon(url('http://mskilab.com/gTrack/hg19/gencode.v19.annotation.gtf.gr.rds')))
if (is.character(gencode))
{
message('pulling GENCODE gtf from ', gencode)
gencode = rtracklayer::import(gencode)
}
cat('loaded gencode rds\n')
tx = gencode[gencode$type =='transcript']
genes = gencode[gencode$type =='gene']
exons = gencode[gencode$type == 'exon']
utr = gencode[gencode$type == 'UTR']
## ut = unlist(utr$tag)
## utix = rep(1:length(utr), sapply(utr$tag, length))
## utr5 = utr[unique(utix[grep('5_UTR',ut)])]
## utr3 = utr[unique(utix[grep('3_UTR',ut)])]
## utr5$type = 'UTR5'
## utr3$type = 'UTR3'
startcodon = gencode[gencode$type == 'start_codon']
stopcodon = gencode[gencode$type == 'stop_codon']
tmp = c(genes, tx, exons, utr, startcodon, stopcodon)
OUT.COLS = c('gene_name', 'gene_id', 'transcript_name', 'transcript_id', 'type', 'exon_number', intersect(names(mcols(tmp)), gencode.cols))
tmp = tmp[, OUT.COLS]
cat('extracted intervals\n')
## compute tx ord of intervals
ord.ix = order(tmp$transcript_id, match(tmp$type, c('gene', 'transcript', 'exon', 'UTR', 'start_codon','stop_codon')))
tmp.rle = rle(tmp$transcript_id[ord.ix])
tmp$tx.ord[ord.ix] = unlist(lapply(tmp.rle$lengths, function(x) 1:x))
tmp = tmp[order(match(tmp$type, c('gene', 'transcript', 'exon', 'UTR', 'start_codon','stop_codon')))]
cat('reordered intervals\n')
tmp.g = tmp[tmp$type != 'transcript']
if (isTRUE(use.gene.ids)){
gencode.composite = split(tmp.g, tmp.g$gene_id)
} else {
gencode.composite = split(tmp.g, tmp.g$gene_name)
}
values(gencode.composite)$id = grl.eval(gencode.composite, gene_name[1])
values(gencode.composite)$gene_sym = grl.eval(gencode.composite, gene_name[1])
if ('gene_label' %in% names(mcols(tmp.g))){
values(gencode.composite)$gene_label = grl.eval(gencode.composite, gene_label[1])
}
if ('gene_label_color' %in% names(mcols(tmp.g))){
values(gencode.composite)$label.color = grl.eval(gencode.composite, gene_label_color[1])
}
if (!grepl("^http",cached.path.collapsed)) {
cat('saving gene collapsed track\n')
saveRDS(gencode.composite, cached.path.collapsed)
}
tmp.t = tmp[tmp$type != 'gene']
gencode.composite = split(tmp.t, tmp.t$transcript_id)
gn = grl.eval(gencode.composite, gene_name[1])
tn = grl.eval(gencode.composite, transcript_name[1])
values(gencode.composite)$id = paste(gn, tn, sep = '-')
values(gencode.composite)$gene_sym = gn
cat('saving transcript track\n')
if (!grepl("^http",cached.path)) {
cat('saving gene collapsed track\n')
saveRDS(gencode.composite, cached.path)
}
genes = NULL
cat(sprintf('cached composite tracks at %s and %s\n', cached.path, cached.path.collapsed))
}
gencode.composite = read.rds.url(ifelse(gene.collapse, cached.path.collapsed, cached.path))
if (drop.rp11)
gencode.composite = gencode.composite[!grepl('^RP11', values(gencode.composite)$gene_sym)]
if (!is.null(genes))
gencode.composite = gencode.composite[values(gencode.composite)$gene_sym %in% genes]
if (!is.null(grep))
{
ix = rep(FALSE, length(gencode.composite))
for (g in grep)
ix = ix | grepl(g, values(gencode.composite)$gene_sym)
gencode.composite = gencode.composite[ix]
}
if (!is.null(grepe))
{
ix = rep(TRUE, length(gencode.composite))
for (g in grepe)
ix = ix & !grepl(g, values(gencode.composite)$gene_sym)
gencode.composite = gencode.composite[ix]
}
cmap = list(type = c(gene = bg.col, transcript = bg.col, exon = cds.col, start_codon = st.col, stop_codon = en.col, UTR = utr.col))
return(suppressWarnings(gTrack(gencode.composite, col = NA, grl.labelfield = 'id', gr.labelfield = 'exon_number',
gr.srt.label = gr.srt.label, cex.label = cex.label,
gr.cex.label = gr.cex.label,
labels.suppress.gr = labels.suppress.gr,
stack.gap = stack.gap, colormaps = cmap, ...)))
}
# @name track.splice
# @title track.splice
#
# Given set of exons and rna bam (eg from tophat) determines junction and exon read density and returns a gTrack object
# of splicing graph
#
# @param ex GRanges of candidate exons
# @param bam path to indexed RNA seq bam
# @param verbose
# @import Rsamtools
# @export
# @author Marcin Imielinski
# track.splice = function(ex = NULL, region = NULL, bam, verbose = TRUE,
# infer.exons = FALSE,
# min.reads = 0, ## only relevant if ex is null or infer.exons = TRUE, here exons are inferred from "N" intervals
# min.exon.width = 10,
# max.exon.width = 1000 ## only relevant if ex is null
# )
# {
# if (!is.null(ex))
# {
# ex = gr.stripstrand(ex)
# ex.ov = gr.findoverlaps(ex, ex)
# }
#
# if (is.null(region))
# reads = read.bam(bam, intervals = ex, pairs.grl = FALSE, verbose = verbose)
# else
# reads = read.bam(bam, intervals = region, pairs.grl = FALSE, verbose = verbose)
#
# if (verbose)
# cat(length(reads), 'reads\n')
#
# if (length(reads)==0)
# {
# if (!is.null(ex))
# {
# ex$expr = 0
# ex$log.expr = log(ex$expr)
# return(gTrack(ex, y.field = 'log.expr', name = 'Log Read Density', height = 30))
# }
# else
# return(gTrack(region[c()]))
# }
#
# sp.reads = splice.cigar(reads, return.grl = FALSE)
#
# if (verbose)
# cat(length(sp.reads), 'spliced fragments\n')
#
# if (is.null(ex))
# infer.exons = TRUE
#
# if (infer.exons)
# {
# sp.reads = sp.reads[seqnames(sp.reads)==seqnames(region) & ranges(sp.reads) %over% ranges(region), ]
# tmp = grdt(sp.reads)
# m.reads = tmp[type != 'N']
# n.reads = tmp[type == 'N']
# ustart = c(min(start(sp.reads)), n.reads[, length(seqnames), by = end][V1>=min.reads, end]+1)
# uend = c(n.reads[, length(seqnames), by = start][V1>=min.reads, start]-1, max(end(sp.reads))) ## 1 after N ends are starts of exons
#
# ustart.maxwidth = m.reads[, max(end-start), keyby = start]
# uend.maxwidth = m.reads[, max(end-start), keyby = end]
#
# new.ex = data.table(seqnames = n.reads$seqnames[1], start = rep(ustart, length(uend)), end = rep(uend, each = length(ustart)))[ (end-start) <= max.exon.width & (end-start) >= min.exon.width, ]
#
# ## we will get exon overload unless we cull a bit
# ## to be parsimonious we only keep enough exons for all the matching parts of reads to "land on"
# ## ie if they start in an exon they will end up in the same one
# ## rather than comp
## uting all the overlaps necessary for this
# ## we just approximate by removing exons for which a smaller one exists that accomodates all of the
# ## reads that start on its start or end on its end
# new.ex[ , maxwidth.start := (end-start) > ustart.maxwidth[list(new.ex$start), V1]]
# new.ex[ , maxwidth.end := (end-start) > uend.maxwidth[list(new.ex$end), V1]]
#
# new.ex[ , shorter.start.exists := !(1:length(seqnames) %in% which.min(end-start)), by = start]
# new.ex[ , shorter.end.exists := !(1:length(seqnames) %in% which.min(end-start)), by = end]
#
# new.ex = seg2gr(new.ex[!shorter.start.exists | !shorter.end.exists, ], seqlengths = seqlengths(sp.reads))[, c()]
#
# new.ex$type = 'inferred'
# if (verbose)
# cat(sprintf('inferred %s total unique exons with min.reads %s, min.exon.width %s, and max.exon.width %s\n', length(new.ex), min.reads, min.exon.width, max.exon.width))
# if (!is.null(ex))
# {
# old.ex = ex;
# ex = sort(unique(gUtils::grbind(ex, new.ex)))
# if (verbose)
# cat(sprintf('Added %s additional exons to yield %s total\n', length(ex)-length(old.ex), length(ex)))
# }
# else
# ex = sort(new.ex)
#
# ex.ov = gr.findoverlaps(ex, ex)
# }
#
# sp.reads = sort(sp.reads[sp.reads$type != 'N'])
# sp.reads = sp.reads[order(sp.reads$rid), ]
#
# tmp = grdt(sp.reads)
# tmp[, spid := 1:length(sp.reads)]
# tmp[, riid := 1:length(seqnames), by = rid]
# setkey(tmp, spid)
# sp.reads$riid = tmp[list(1:length(sp.reads)), riid]
#
# ov = gr.findoverlaps(ex, sp.reads, scol = c('rid', 'riid'), verbose = verbose, return.type = 'data.table')
# if (length(ov)==0)
# {
# ex$expr = 0
# ex$log.expr = log(ex$expr)
# return(gTrack(ex, y.field = 'log.expr', name = 'Log Read Density', height = 30))
# }
# ov[, width := end-start]
# ov[, match.left.exon := as.numeric(start == start(ex)[query.id])]
# ov[, match.right.exon := as.numeric(end == end(ex)[query.id])]
# ov[, match.left.read := as.numeric(start == start(sp.reads)[subject.id])]
# ov[, match.right.read := as.numeric(end == end(sp.reads)[subject.id])]
# ov = ov[(match.left.read | match.left.exon) & (match.right.read | match.right.exon), ]
# setkeyv(ov, c('rid', 'riid'))
#
# sp.rid = unique(ov$rid[ov$riid>1])
# ij = ov[list(sp.rid), list(from = rep(query.id, each = length(query.id)), to = rep(query.id, length(query.id)),
# val = width[rep(1:length(query.id), each = length(query.id))] + width[rep(1:length(query.id), length(query.id))],
# from.riid = rep(riid*match.right.exon*match.right.read, each = length(query.id)),
# to.riid = rep(riid*match.left.exon*match.left.read, length(query.id))), by = rid]
#
# ij = ij[(ij$to.riid - ij$from.riid) == 1 & ij$to.riid !=0 & ij$from.riid != 0, ]
# setkeyv(ij, c('from', 'to'))
# ij = ij[!list(ex.ov$query.id, ex.ov$subject.id), ]
#
# # edges = ij[, list(val = sum(val)/(width(ex)[from] + width(ex)[to])), keyby = list(from, to)]
# edges = ij[, list(val = sum(val)), keyby = list(from, to)]
# #edges = edges[, val := 0]
# edges = edges[, lwd := affine.map(log(val+1), c(0, 6), cap = TRUE)]
# edges = edges[, v := 5]
# edges = edges[, h := 2]
# edges = edges[, col := alpha('gray10', affine.map(log(val+1), c(0.1,0.8), cap = TRUE))]
# edges = edges[, cex.arrow := 0]
# ex$expr = ov[(match.right.exon*match.right.read + match.left.exon*match.right.read) | (match.left.read==1 & match.right.read==1), sum(width), keyby = query.id][list(1:length(ex)), V1]/width(ex)
#
# ex$log.expr = round(log10(ex$expr), 1)
# ex$ywid = 0.8
# ex$border = 'black'
# ex$col = alpha('blue', 0.4)
# td.ex = gTrack(ex, y.field = 'log.expr', edges = edges, name = 'Log Read Density', height = 30)
# return(td.ex)
# }
##########################
#' @name draw.ranges
#' @title draw.ranges
#' @description
#'
#' Draws (genomic, Iranges, Rle, or df with start and end field) intervals
#' as rectangles of thickness "lwd", col "col".
#'
#' If "col" or "y" are set to null then will examine "x" to see if the respective fields are contained within (e.g. as RangedData or data.frame)
#'
#' Adds optional backbone between beginning of first range
#' and end of last on existing plot
#' @keywords internal
#' @author Marcin Imielinski
draw.ranges = function(x, y = NULL, lwd = 0.5, col = "black", border = col, label = NULL,
draw.backbone = F, # if TRUE lines will be drawn linking all ranges at a single y-level +/- (a single y x group level if group is specified),
group = NULL, # a vector of group labels same size as x (char or numeric), non-NA groups with more than 1 member will be connected by backbone (if draw.backbone = T)
col.backbone = col[1],
cex.label = 1,
srt.label = 45,
adj.label = c(0.5, 0),
angle, # vector of angles (0 = rectangle -45 leftward barb, +45 rightward barb
circles = FALSE, # if TRUE will draw circles at range midpoints instead of polygons
bars = FALSE, # if TRUE will draw bars from y0.bar to the y locations
points = NA, # if integer then will draw points at the midpoint with pch value "points"
lines = FALSE, # if TRUE will connect endpoints of ranges with lines
y0.bar = NULL, # only relevant if bars = T, will default to pars('usr')[3] if not specified
strand.dir = T, # if so, will interpret strand field as "direction sign" for interval, + -> right, - -> left, only makes difference if draw.pencils = T
xlim = NULL, # only relevant if new.plot = T, can be ranges object
ylim = NULL, # only relevant if new.plot = T
clip = NULL, # GRanges / IRanges or 1 row df with $start, $end denoting single region to clip to .. will not draw outside of this region
lwd.border = 1,
lwd.backbone = 1,
verbose=FALSE,
...)
{
PAD = 0 ##0.5;
if (is.null(y))
#if (is.data.frame(x))
# y = stack_segs(x)
#else
y = IRanges::disjointBins(IRanges::ranges(x))
#if (inherits(x, 'Rle'))
# x = RangedData(x)
if (inherits(x, 'GappedAlignments'))
x = as(x, 'GRanges')
if (inherits(x, 'Ranges') | inherits(x, 'GRanges') | inherits(x, 'RangedData'))
x = standardize_segs(x)
if (nrow(x)==0)
return()
if (!is.null(y))
x$y = y;
x$group = group;
if (!is.null(col) & is.null(x$col))
x$col = col;
if (!is.null(border) & is.null(x$border))
x$border = border;
if (!is.null(lwd.border) & is.null(x$lwd.border))
x$lwd.border = lwd.border;
if (!is.null(label) & is.null(x$label))
x$label = label;
if (!is.null(circles) & is.null(x$circles))
x$circles = circles
if (!is.null(bars) & is.null(x$bars))
x$bars = bars
if (is.null(x$points))
x$points = NA
if (!is.na(points) & all(is.na(x$points)))
x$points = points
if (is.null(x$lines))
x$lines = NA
if (!is.na(lines) & all(is.na(x$lines)))
x$lines = lines
if (!is.null(clip))
{
clip = standardize_segs(clip);
x = clip.seg(x, clip)
}
if (is.null(srt.label))
srt.label = NA
if (is.na(srt.label))
srt.label = 45
if (strand.dir & !is.null(strand))
x$direction = x$strand
x$lwd = lwd
if (nrow(x)>0)
{
x$group = paste(x$group, x$y);
if (draw.backbone)
{
gcount = table(x$group);
x2 = x[x$group %in% names(gcount)[gcount>1], ]
if (nrow(x2)>0)
{
b.st = aggregate(pos1 ~ group, data = x2, FUN = min)
b.en = aggregate(pos2 ~ group, data = x2, FUN = max)
b.y = aggregate(y ~ group, data = x2, FUN = function(x) x[1])
seg.coord = data.frame(x0 = b.st[,2], x1 = b.en[,2], y = b.y[,2]);
if (!is.null(clip))
{
seg.coord$x0 = pmax(seg.coord$x0, clip$pos1)
seg.coord$x1 = pmin(seg.coord$x1, clip$pos2)
}
if (verbose)
print('Finished computing backbones')
segments(seg.coord$x0, seg.coord$y, seg.coord$x1, seg.coord$y, col = col.backbone, lwd = lwd.backbone);
}
}
if (is.null(x$direction))
x$direction = "*"
## if (is.na(circles))
## circles = F
## if (is.na(lines))
## lines = F
if (is.null(x$circles))
x$circles = FALSE
if (is.null(x$points))
x$points = FALSE
if (is.null(x$lines))
x$lines = FALSE
if (is.null(x$bars))
x$bars = FALSE
x$bars = ifelse(is.na(x$bars), FALSE, x$bars)
x$circles = ifelse(is.na(x$circles), FALSE, x$circles)
x$lines = ifelse(is.na(x$lines), FALSE, x$lines)
x$points = ifelse(is.na(x$points), FALSE, x$points)
x$rect = !x$points & !x$bars & !x$circles & !x$lines ## default barbs are only drawn if points / circles / lines/ bars not specified
if (any(x$rect))
{
x.tmp = x[which(x$rect), ]
main.args = list(x0 = x.tmp$pos1 - PAD, y0 = x.tmp$y-x.tmp$lwd/2, x1 = x.tmp$pos2 + PAD, y1 = x.tmp$y + x.tmp$lwd/2, col = x.tmp$col, border = x.tmp$border , angle = angle*c('*'=0, '+'=1, '-'=-1)[x.tmp$direction], lwd = x.tmp$lwd.border)
other.args = list(...);
other.args = other.args[setdiff(names(other.args), names(main.args))]
do.call(barbs, c(main.args, other.args))
}
if (any(x$points))
{
x.tmp = x[which(x$points), ]
points((x.tmp$pos1 + x.tmp$pos2)/2, x.tmp$y, pch = points, border = x.tmp$border, col = x.tmp$col, cex = x.tmp$lwd.border)
# points((x$pos1 + x$pos2)/2, x$y-x$lwd/2, pch = points, border = x$border, col = x$col, cex = x$lwd.border)
}
if (any(x$circles))
{
x.tmp = x[which(x$circles), ]
points((x.tmp$pos1 + x.tmp$pos2)/2, x.tmp$y, pch = 19, col = x.tmp$col, cex = x.tmp$lwd.border)
points((x.tmp$pos1 + x.tmp$pos2)/2, x.tmp$y, pch = 1, col = x.tmp$border, cex = x.tmp$lwd.border)
## points((x$pos1 + x$pos2)/2, x$y, pch = points, col = x$col, cex = x$lwd.border)
## points((x$pos1 + x$pos2)/2, x$y-x$lwd/2, pch = points, col = x$col, cex = x$lwd.border)
## points((x$pos1 + x$pos2)/2, x$y-x$lwd/2, pch = 1, col = x$border, cex = x$lwd.border)
}
if (any(x$bars))
{
if (is.null(y0.bar))
y0.bar = par('usr')[3]
x.tmp = x[which(x$bars), ]
rect(x.tmp$pos1, ifelse(y0.bar<x.tmp$y, rep(y0.bar, nrow(x)), x.tmp$y), x.tmp$pos2, ifelse(y0.bar<x.tmp$y, x.tmp$y, rep(y0.bar, nrow(x))),
col = x.tmp$col, border = x.tmp$col, lwd = x.tmp$lwd.border/2)
# rect(x$pos1, rep(y0.bar, nrow(x)), x$pos2, x$y, col = x$col, border = NA, lwd = NA)
}
if (any(x$lines))
{
x.tmp = x[which(x$lines), ]
lines(as.vector(t(cbind(x.tmp$pos1, x.tmp$pos2))), as.vector(t(cbind(x.tmp$y, x.tmp$y))), col = as.vector(t(cbind(x.tmp$col, x.tmp$col))),
lwd = as.vector(t(cbind(x.tmp$lwd.border, x.tmp$lwd.border))))
}
## if (any(!is.na(x$points)))
## {
## points((x$pos1 + x$pos2)/2, x$y-x$lwd/2, pch = x$points, col = x$col, cex = x$lwd.border)
## }
## if (circles)
## {
## points((x$pos1 + x$pos2)/2, x$y-x$lwd/2, pch = 19, col = x$col, cex = x$lwd.border)
## points((x$pos1 + x$pos2)/2, x$y-x$lwd/2, pch = 1, col = x$border, cex = x$lwd.border)
## }
if (!is.null(x$label))
{
has.label = !is.na(x$label)
if (any(has.label))
{
if (nrow(x)>1)
OFFSET = min(diff(sort(x$y)))/2
else
OFFSET = diff(par('usr')[3:4])/100
if (adj.label[2] == 1)
{
text(rowMeans(x[has.label, c('pos1', 'pos2')]), x$y[has.label]+x$lwd[has.label]/2+OFFSET, as.character(x$label[has.label]),
adj = adj.label, cex = 0.5*cex.label, srt = srt.label)
}
else if (adj.label[2] == 0.5)
{
text(rowMeans(x[has.label, c('pos1', 'pos2')]), x$y[has.label], as.character(x$label[has.label]),
adj = adj.label, cex = 0.5*cex.label, srt = srt.label)
}
else
{
text(rowMeans(x[has.label, c('pos1', 'pos2')]), x$y[has.label]-x$lwd[has.label]/2-OFFSET, as.character(x$label[has.label]),
adj = adj.label, cex = 0.5*cex.label, srt = srt.label)
}
}
}
}
}
#' @name draw.grl
#' @title draw.grl
#' @description
#'
#' wrapper function around draw.ranges to draw GRangesList across a set of genomic windows
#' By default, these windows correspond to the smallest set of "covered" regions (ie containing at least one
#' grl in their span, but can also be specified as an argument.
#'
#' basically same idea as draw.granges, except can more compactly draw / stack paired / group reads or
#' gene models / alignments and deal with "grouped" intervals
#'
#' additional values features of input grl that can impact drawing:
#' ywid = vertical thickness of y segments
#' border.col = color border (can override colormap)
#'
#' if optional arg draw.var == T
#' will either call varbase on grl (e.g. if grl is a GRanges or GappedAlignment object)
#' to get variant ranges or will take vars object (GRangesList that is of same length as grl)
#' and will display variant ranges according to color scheme in arg "var.col"
#' "var.col" is a named vector that maps variant types to colors, with the following syntax:
#'
#' XA, XT, XG, XC --> colors of base substitutions
#' S, D, I --> colors of sooutft clipped, deletions, insertions
#' e.g.
#' var.col = c(XA = 'blue', S = 'yellow', 'I' = green') would replace the default variant coloring of
#' "A" substitutions, soft clipped regions, and insertions with blue, yellow, and green colors, respectively.
#'
#' @keywords internal
#' @author Marcin IMielinski
#' @importFrom GenomicRanges GRanges values ranges width strand values<- strand<- seqnames coverage ranges<-
#' @importFrom data.table := setkeyv
#' @importFrom GenomeInfoDb Seqinfo seqinfo keepSeqlevels seqlevels seqlengths seqlevels<- seqlengths<- genome<- seqnames
draw.grl = function(grl,
y = NULL, # can be either vector of length grl, or data.frame row / list with fields $start and $end
# specifying y coordinates to "pack" the granges into (or just length 2 list)
# note this is different from ylim, which determines the size of the canvas
ylim = NULL, # if null and y provided, will find range of y and plot
ylim.subplot = NULL,
ywid = NULL,
edges = NULL, ## data.frame specifying edges connecting items of grl, with fields $from $to indexing grl and optional fields $lwd, $col, $lty specifying formatting options for connectors, for gr the from arrow will leave the right side of a + range and left side of a - range, and enter the left side of a + range and right side of a - range. For grl, from arrows will be drawn from the last range in the "from" list item to the first range in the "to" list item
draw.paths = F, # if draw.paths = T will treat grl's as paths / contigs,
# connecting intervals joined by arrowed curved arcs using alternate stacking algo (ignoring y information)
draw.var = F, # if true, then varbase will be called on grl or unlist(grl)
var = NULL, # optional GRangesList of same length as grl specifying variant bases with value cols $type, $varbase (ie output of varbase)
# corresponding to each grl item
var.col = NULL, # named vector with variant colors to override - valid names are XA, XT, XG, XC, S, D, I
var.soft = T, ## whether to draw soft clips for varbase
windows, # windows specifies windows, can have optional meta data features $col and $border
win.gap = NULL,
stack.gap,
min.gapwidth = 1, # only applies if windows is not specified
col = NULL, border = NA, # all of these variables can be scalar or vectors of length(grl),
# can be over-ridden by values / columns in grl
col.backbone = alpha('gray', 0.8),
gr.colorfield = NULL, ## values field in the gr from which colors can be mapped
gr.colormap = NULL, ## named vector mapping fields in the gr.colorfield to colors, if unspecified brewer.master() will be applied
gr.labelfield = NULL, ## field of gr labels to draw.
grl.labelfield = NULL, ## field of grl to draw as label
grl.colorfield = NULL, ## field of grl to set colors for labels
grl.cexfield = NULL, ## field of grl to set cex for labels
leg.params,
labels = NULL, # vector of length(grl)
labels.suppress = F,
labels.suppress.grl = labels.suppress,
labels.suppress.gr = labels.suppress,
spc.label = 0.05, # number between 0 and 1 indicating spacing of label
adj.label = c(0, 1),
cex.label = 1,
gr.cex.label = 0.8 * cex.label,
gr.srt.label = 0,
gr.adj.label = c(0,0.5),
new.plot, new.axis,
sep.lty = 2,
sep.lwd = 1,
sep.bg.col = 'gray95',
sep.draw = TRUE,
y.pad, # this is the fractional padding to put on top and bottom of ranges if y is specified as $start and $end pair (def was 0.05)
xaxis.prefix = '', xaxis.suffix = 'MB', xaxis.unit = 1, xaxis.round = 3,
xaxis.interval = 'auto', xaxis.pos = 1,
xaxis.pos.label, xaxis.cex.label,
xaxis.newline = FALSE,
xaxis.chronly = FALSE,
xaxis.ticklen = 1,
xaxis.width = TRUE,
xaxis.label.angle = 0,
xaxis.cex.tick = 1,
ylim.grid = ylim, # can be also specified directly for plots with multiple axes and/or non-numeric tracks
y.grid = NA, # if non NA, then the number specifies the spacing between y grid (drawn from ylim[1] to ylim[2]), can also be named vector mapping grid.tick labels to plot locations
ylab = NULL,
y.grid.col = alpha('gray', 0.5),
y.grid.pretty = 5,
y.grid.cex = 1,
y.grid.lty = 2,
y.grid.lwd = 1,
path.col = 'black',
path.lwd = 1,
path.col.arrow = path.col,
path.cex.arrow = 1,
path.stack.y.gap = 1,
path.stack.x.gap = 0,
path.cex.v = 1,
path.cex.h = 1,
draw.backbone = NULL,
xlim = c(0, 20000), # xlim of canvas
points = NA, ## if non NA then will draw a given point with pch style
circles = FALSE, ## only one of these should be true, however if multiple are true then they will be interpreted in this order
bars = FALSE,
y0.bar = NULL,
lines = F,
angle, # angle of barbs to indicate directionality of ranges
verbose=FALSE,
triangle=FALSE, # trigger a triangle matrix plot
ylim.parent=NULL, ## ylim of the full thing. This is importat for angle preseveration
legend.params = list(plot=TRUE),
bg.col = NA, ## background of whole plot
...)
{
now = Sys.time();
## PATCH: we are forgetting about any ylim.subplot settings done above .. WHY?
# ylim.subplot = NULL
empty.plot = FALSE
## PATCH: last minute defaults
if (is.na(xaxis.width))
xaxis.width = TRUE
if (is.na(xaxis.chronly))
xaxis.chronly = FALSE
if (is.na(xaxis.label.angle))
xaxis.label.angle = 0
if (length(grl)>0)
{
if (is.null(draw.backbone))
draw.backbone = TRUE
# if (inherits(grl, 'GappedAlignments'))
# grl = ga2gr(grl)
if ((inherits(grl, 'GRanges')))
{
if (!is.null(windows)) ## hack to get over stupid GRanges speed issues when we have a giant GRanges input (eg coverage)
{
strand(windows) <- rep("*", length(windows)) ## don't need strand on windows, mess up intersections
ix <- gUtils::gr.in(grl, windows)
grl = grl[ix]
if (!is.null(col))
if (length(col)!=1)
col = col[ix]
if (!is.null(border))
if (length(border)!=1)
border = border[ix]
if (!is.null(ywid))
if (length(ywid)!=1)
ywid = ywid[ix]
if (!is.null(y))
if (!is.list(y))
if (length(y)!=1)
y = y[ix]
if (!is.null(labels))
if (!is.list(labels))
if (length(labels)!=1)
labels = labels[ix]
if (!is.null(edges))
if (nrow(edges)>0 & all(c('from', 'to') %in% colnames(edges)))
{
if (data.table::is.data.table(edges))
edges = as.data.frame(edges)
ix.i = which(ix)
edges = edges[edges$from %in% ix.i | edges$to %in% ix.i,];
edges$from = match(edges$from, ix.i)
edges$to = match(edges$to, ix.i)
}
}
gr.names = names(grl);
names(grl) = NULL;
if (length(grl)>0)
grl = GenomicRanges::split(grl, 1:length(grl))
else
grl = GenomicRanges::GRangesList(grl)
names(grl) = gr.names;
}
if (is.null(labels))
if (is.null(values(grl)$labels))
labels = names(grl)
else
labels = values(grl)$labels
grl.colors = NA
if (is.null(grl.colorfield) || is.na(grl.colorfield))
if (!is.null(values(grl)$label.color))
grl.colors = values(grl)$label.color
grl.cex = cex.label
if (is.null(grl.cexfield) || is.na(grl.cexfield))
if (!is.null(values(grl)$cex.label))
grl.cex = values(grl)$cex.label
# replace missing or non-numeric values with cex.label
grl.cex[which(is.null(grl.cex) | is.na(grl.cex) | !is.numeric(grl.cex))] = cex.label
# make sure names are unique
names(grl) = 1:length(grl);
if (is.na(labels.suppress.grl))
labels.suppress.grl = FALSE
if (is.na(labels.suppress.gr))
labels.suppress.gr = FALSE
if (is.na(draw.var))
draw.var = F
if (is.na(draw.paths))
draw.paths = F
if (!is.null(gr.colorfield))
if (all(is.na(gr.colorfield)))
gr.colorfield = NULL
if (!is.null(gr.colormap))
if (all(is.na(gr.colormap)))
gr.colormap = NULL
else if (is.null(names(gr.colormap)))
gr.colormap = NULL
if (!is.null(col))
if (all(is.na(col)))
col = NULL
## postpone further color processing until later
if (is.null(col))
if (!is.null(values(grl)$col))
col = values(grl)$col
else if (is.null(gr.colorfield) & is.null(gr.colormap))
col = NA # 'black'
else
col = NA
if (is.na(col.backbone))
col.backbone = NULL
if (is.null(col.backbone))
col.backbone = col
# if (is.na(border))
# border = NULL
if (is.null(border))
border = NA
grl.props = cbind(data.frame(group = names(grl), col = col, stringsAsFactors = F), as.data.frame(values(grl)))
grl.props$border = border;
grl.props$ywid = ywid;
if (!is.null(y) & !is.list(y) & length(y)>0) ## if y coordinate is specified for the ranges
{
grl.props$y = y;
bad = is.na(y)
bad[which(is.infinite(y))] = TRUE
if (any(bad))
{
grl.props = grl.props[!bad, ]
grl = grl[!bad]
if (!is.null(labels))
labels = labels[!bad]
y = grl.props$y
}
}
if (!is.null(grl.labelfield))
{
if (!is.na(grl.labelfield))
{
if (grl.labelfield %in% names(grl.props))
grl.props$grl.labels = grl.props[, grl.labelfield]
}
else if (!is.null(labels))
grl.props$grl.labels = labels ## use $grl.labels to allow labeling of individual grs
}
else if (!is.null(labels))
if (!is.na(labels[1])) ## will only be null if labels is NULL and names(grl) was NULL
grl.props$grl.labels = labels ## use $grl.labels to allow labeling of individual grs
grl.props$grl.cols = grl.colors
if (!is.null(grl.props$grl.labels) && !is.null(grl.colorfield) && grl.colorfield %in% names(grl.props)){
grl.props$grl.cols = grl.props[, grl.colorfield] ## use $grl.cols to allow custom coloring of individual labels
}
grl.props$grl.cex = grl.cex
if (!is.null(grl.props$grl.labels) && !is.null(grl.cexfield) && grl.cexfield %in% names(grl.props)){
grl.props$grl.cex = grl.props[, grl.cexfield] ## use $grl.cex to allow custom sizing of individual labels
}
gr = tryCatch(grl.unlist(grl), error = function(e)
{
## ghetto solution if we get GRanges names snafu
gr = unlist(grl);
if (length(gr)>0)
{
tmpc = textConnection(names(gr));
cat('budget .. \n')
gr$grl.ix = read.delim(tmpc, sep = '.', header = F)[,1];
# gr$grl.iix = levapply(rep(1, length(gr)), gr$grl.ix, FUN = function(x) 1:length(x))
gr$grl.iix = data.table::data.table(ix = gr$grl.ix)[, iix := 1:length(ix), by = ix][, iix]
close(tmpc)
}
return(gr)
})
gr$group = grl.props$group[gr$grl.ix]
gr$group.ord = gr$grl.iix
gr$first = gr$grl.iix == 1
last = iix = NULL ## NOTE fix
if (length(gr)>0)
gr$last = data.table::data.table(iix = as.numeric(gr$grl.iix), ix = gr$grl.ix)[, last := iix == max(iix), by = ix][, last]
# gr$last = levapply(gr$grl.iix, gr$grl.ix, FUN = function(x) x == max(x))
grl.props$group = as.character(grl.props$group)
S4Vectors::values(gr) = cbind(as.data.frame(values(gr)),
grl.props[match(values(gr)$group, grl.props$group), setdiff(colnames(grl.props), c(colnames(values(gr)), 'group', 'labels')), drop = FALSE])
#####
## variant drawing
####
if (draw.var & is.null(var) & length(gr)>0)
{
# var = bamUtils::varbase(gr[gr.in(gr, windows)], soft = var.soft)
gix = which(gr.in(gr, windows))
var = varbase(gr[gix], soft = var.soft)
if (any(iix <- var$type == 'I'))
end(var[ix]) = end(var[ix])+1
names(var) = gix
}
else
gix = NULL
if (!is.null(var))
if (inherits(var, 'GRangesList'))
{
VAR.COL = c('XA' = 'green', 'XG' = 'brown', 'XC' = 'blue', 'XT' = 'red', 'D' = 'white', 'I'= 'purple', 'N' = alpha('gray', 0.2), 'XX' = 'black', 'S' = alpha('pink', 0.9))
if (!is.null(var.col))
VAR.COL[names(var.col)] = var.col;
names(var) = NULL
if (!is.null(gix))
names(var) = gix
else
names(var) = 1:length(var)
# var.group = as.numeric(as.data.frame(var)$element)
var.gr = grl.unlist(var)
if (length(var.gr)>0)
{
var.group = names(var)[var.gr$grl.ix]
# var.gr = gr.stripstrand(unlist(var))
# inherit properties from gr
values(var.gr) = cbind(as.data.frame(values(var.gr)),
as.data.frame(values(gr)[as.numeric(var.group), setdiff(names(values(gr)), c('labels'))]))
if (!is.null(gr.labelfield))
if (!is.na(gr.labelfield))
values(var.gr)[, gr.labelfield] = NA
var.gr$col.sig = as.character(var.gr$type);
xix = var.gr$type == 'X'
var.gr$col.sig[xix] = paste(var.gr$col.sig[xix], var.gr$varbase[xix], sep = '')
var.gr$col = VAR.COL[var.gr$col.sig]
var.gr$border = var.gr$col
var.gr$first = FALSE
var.gr$last = FALSE
if (draw.paths) ## if we are drawing paths, then need to setdiff variant vs non variant parts of edges and re-order
{
## remove soft clips
var.gr = var.gr[!(var.gr$type %in% c('H', 'S'))]
gr2 = gr;
GenomicRanges::strand(var.gr) == GenomicRanges::strand(gr)[var.gr$grl.ix]
gr$grl.iix = as.numeric(gr$grl.iix)
## now doing some ranges acrobatics to find all the pieces of gr that are not in var.gr
## TODO: speed up, remove awkawardness
ir = IRanges::ranges(gr)
var.ir = IRanges::ranges(var.gr)
tmp.ix = split(1:length(var.gr), var.gr$grl.ix)
tmp.l = lapply(names(tmp.ix), function(i) IRanges::disjoin(c(ir[as.numeric(i)], var.ir[tmp.ix[[i]]])))
tmp.ogix = rep(as.numeric(names(tmp.ix)), sapply(tmp.l, length))
tmp.ir = do.call(c, tmp.l)
tmp.gr = GenomicRanges::GRanges(seqnames(gr)[tmp.ogix], tmp.ir, seqlengths = GenomeInfoDb::seqlengths(gr), og.ix = tmp.ogix)
tmp.ov = gr.findoverlaps(tmp.gr, var.gr)
tmp.ov = tmp.ov[tmp.gr$og.ix[tmp.ov$query.id] == var.gr$grl.ix[tmp.ov$subject.id] ]
new.gr = tmp.gr[!(1:length(tmp.gr) %in% tmp.ov$query.id), ] ## only keep the non variant pieces
GenomicRanges::strand(new.gr) = GenomicRanges::strand(gr)[new.gr$og.ix]
values(new.gr) = cbind(as.data.frame(values(gr)[new.gr$og.ix, ]) , og.ix = new.gr$og.ix)
var.gr$og.ix = var.gr$grl.ix
GenomicRanges::strand(var.gr) = GenomicRanges::strand(gr)[var.gr$og.ix]
# var.gr$group = as.numeric(as.character(gr$group[var.gr$og.ix]))
new.gr = grbind(new.gr, var.gr, gr[setdiff(1:length(gr), var.gr$grl.ix)])
new.gr$grl.iix = as.numeric(gr$grl.iix[new.gr$og.ix])
new.ord = mapply(function(x, y, z) if (y[1]) x[order(z)] else rev(x[order(z)]),
split(1:length(new.gr), new.gr$og.ix), split(as.logical(GenomicRanges::strand(new.gr)=='+'), new.gr$og.ix), split(start(new.gr), new.gr$og.ix))
new.ix = unlist(lapply(new.ord, function(x) ((1:length(x))-1)/length(x)))
new.gr$grl.iix[unlist(new.ord)] = new.gr$grl.iix[unlist(new.ord)] + new.ix
gr = new.gr[order(new.gr$group, new.gr$grl.iix), ]
}
gr = grbind(gr, var.gr)
}
else
{
gr = grbind(gr, var.gr)
}
}
if (length(gr)>0)
names(gr) = 1:length(gr)
if (is.null(windows)) ## find covered windows in provided grl
{
seqlevels(gr) = seqlevels(gr)[seqlevels(gr) %in% as.character(seqnames(gr))]
windows = as(coverage(gr), 'GRanges');
windows = windows[values(windows)$score!=0]
windows = reduce(windows, min.gapwidth = min.gapwidth);
}
else if (!is(windows, 'GRanges'))
if (is(windows, 'GRangesList'))
windows = unlist(windows)
else ## assume it's a seqinfo object or an object that has a seq
windows = si2gr(windows)
if (is.null(win.gap))
win.gap = mean(width(windows))*0.2
if (sum(as.numeric(width(windows)))==0)
stop('Windows have width 0')
if (verbose) {
print('Before flatmap')
print(Sys.time() - now)
}
## add 1 bp to end for visualization .. ranges avoids weird width < 0 error
if (length(gr)>0)
{
IRanges::ranges(gr) = IRanges::IRanges(start(gr), pmax(end(gr), pmin(end(gr)+1, GenomeInfoDb::seqlengths(gr)[as.character(seqnames(gr))], na.rm = T), na.rm = T)) ## jeremiah commented
# end(gr) = pmax(end(gr), pmin(end(gr)+1, seqlengths(gr)[as.character(seqnames(gr))], na.rm = T), na.rm = T)
}
suppressWarnings(end(windows) <- end(windows) + 1) ## shift one needed bc gr.flatmap has continuous convention, we have categorical (e.g. 2 bases is width 2, not 1)
mapped = gr.flatmap(gr, windows, win.gap);
grl.segs = mapped$grl.segs;
window.segs = mapped$window.segs;
if (verbose) {
print('After flatmap')
print(Sys.time() - now)
}
grl.segs$border = as.character(grl.segs$border)
grl.segs$col = as.character(grl.segs$col)
grl.segs$group = as.character(grl.segs$group)
grl.segs$strand = as.character(grl.segs$strand)
if (!is.null(gr.labelfield))
if (!is.na(gr.labelfield))
if (gr.labelfield %in% names(grl.segs))
grl.segs$label = grl.segs[, gr.labelfield]
if (nrow(grl.segs)==0)
{
# warning('No ranges intersecting window');
# return()
}
if (is.list(y))
{
if (all(c('start', 'end') %in% names(y)))
ylim.subplot = c(y$start[1], y$end[1])
else
ylim.subplot = c(y[[1]], y[[length(y)]])
}
if (nrow(grl.segs)>0)
{
if (!draw.paths)
{
if (is.null(y) | !is.null(ylim.subplot))
{
pos1 = aggregate(pos1 ~ group, data = grl.segs, FUN = min);
pos2 = aggregate(pos2 ~ group, data = grl.segs, FUN = max);
pos1 = structure(pos1[,2], names = pos1[,1]) - round(stack.gap/2);
pos2 = structure(pos2[,2]-1, names = pos2[,1]) + round(stack.gap/2);
ix = order(as.numeric(names(pos1)))
pos1 = pos1[ix]
pos2 = pos2[ix]
## FIX HERE .. avoid using disjointBins
## these gymnastics allow us to use disjointBins (which works on IRanges)
## without getting integer overflow
if (max(c(pos1, pos2))>2e9)
{
r = range(c(pos1, pos2))
pos1 = affine.map(pos1, xlim = r, ylim = c(0, 2e9))
pos2 = affine.map(pos2, xlim = r, ylim = c(0, 2e9))
## m = max(c(pos1,pos2));
## pos1 = ceiling(pos1/m*2e9)
## pos2 = floor(pos2/m*2e9);
}
# bin ranges
y.bin = IRanges::disjointBins(IRanges::IRanges(pos1, pos2))
m.y.bin = max(y.bin)
if (is.null(ylim))
ylim = c(1, m.y.bin) + c(-0.5*m.y.bin, 0.5*m.y.bin)
## squeeze y coordinates into provided (or inferred) ylim
if (!is.null(ylim.subplot))
tmp.ylim = ylim.subplot
else
tmp.ylim = ylim
## provide bottom and top padding of y.bin
y.pad = max(c(0, min(0.49, y.pad)));
# tmp.ylim = tmp.ylim + c(1, -1)*y.pad*diff(tmp.ylim);
y.pad = pmin(1/(m.y.bin+1)/2, 0.125)
tmp.ylim = tmp.ylim + c(1, -1)*y.pad*diff(tmp.ylim);
y = structure(affine.map(y.bin, tmp.ylim), names = names(pos1));
grl.segs$y = y[grl.segs$group]
}
else ## data is numeric, i.e. has some kind of y data
{
if (is.null(ylim))
if (any(!is.na(y[!is.infinite(y)])))
{
tmp.ylim = range(y[!is.infinite(y)], na.rm = T)
ylim = tmp.ylim + c(-1, 1)*0.2*diff(tmp.ylim) + c(-1, 0.2)*y.pad*diff(tmp.ylim)
}
else
ylim = c(0,10)
}
}
else ## draw.paths = T --> will treat each grl as a sequence, which will be joined by connectors
{
ix.l = lapply(split(1:nrow(grl.segs), grl.segs$group), function(x) x[order(grl.segs$group.ord[x])])
grl.segs$y.relbin = NA
## we want to layout paths so that we prevent collisions between different paths
grl.segs$y.relbin[unlist(ix.l)] = unlist(lapply(ix.l, function(ix)
{
# find runs where start[i+1]>end[i] and strand[i] == strand[i+1] = '+'
# and end[i+1]<start[i] and strand[i] == strand[i+1] = '-'
if (length(ix)>1)
{
iix = 1:(length(ix)-1)
concordant = ((grl.segs$pos1[ix[iix+1]] >= grl.segs$pos2[ix[iix]]
& grl.segs$strand[ix[iix+1]] != '-' & grl.segs$strand[ix[iix]] != '-') |
(grl.segs$pos2[ix[iix+1]] <= grl.segs$pos1[ix[iix]]
& grl.segs$strand[ix[iix+1]] == '-' & grl.segs$strand[ix[iix]] == '-'))
return(c(0, cumsum(!concordant)))
}
else
return(0)
}))
contig.lim = data.frame(
group = names(vaggregate(y.relbin ~ group, data = grl.segs, FUN = max)),
pos1 = vaggregate(pos1 ~ group, data = grl.segs, FUN = min) - round(stack.gap)/2,
pos2 = vaggregate(pos2 ~ group, data = grl.segs, FUN = max) + round(stack.gap)/2,
height = vaggregate(y.relbin ~ group, data = grl.segs, FUN = max)
);
contig.lim$width = contig.lim$pos2 - contig.lim$pos1
contig.lim$y.bin = 0;
contig.lim = contig.lim[order(-contig.lim$width), ]
if (nrow(contig.lim)>1)
for (i in 2:nrow(contig.lim))
{
# find lowest level at which there is no clash with this and previously stacked segments
if (max(c(contig.lim$pos1, contig.lim$pos2))>2e9)
{
m = max(c(contig.lim$pos1,contig.lim$pos2));
contig.lim$pos1 = ceiling(contig.lim$pos1/m*2e9)
contig.lim$pos2 = floor(contig.lim$pos2/m*2e9);
}
ir1 = IRanges::IRanges(contig.lim[1:(i-1), 'pos1'], contig.lim[1:(i-1), 'pos2'])
ir2 = IRanges::IRanges(contig.lim[i, 'pos1'], contig.lim[i, 'pos2'])
# clash = which(gUtils::gr.in(ir1, ir2 + path.stack.x.gap))
# clash = which(gUtils::gr.in(ir1, ir2 + path.stack.x.gap))
clash = which(ir1 %over% (ir2 + path.stack.x.gap))
pick = clash[which.max(contig.lim$y.bin[clash] + contig.lim$height[clash])]
contig.lim$y.bin[i] = c(contig.lim$y.bin[pick] + contig.lim$height[pick] + path.stack.y.gap, 0)[1]
}
grl.segs$y.bin = contig.lim$y.bin[match(grl.segs$group, contig.lim$group)] + grl.segs$y.relbin + 1
m.y.bin = max(grl.segs$y.bin)
if (is.null(ylim))
ylim = c(1, m.y.bin) + c(-0.5*m.y.bin, 0.5*m.y.bin)
## squeeze y coordinates into provided (or inferred) ylim
if (!is.null(ylim.subplot))
tmp.ylim = ylim.subplot
else
tmp.ylim = ylim
## provide bottom and top padding of y.bin
#
y.pad = 1/(m.y.bin+1)/2
y.pad = pmin(1/(m.y.bin+1)/2, 0.125)
tmp.ylim = tmp.ylim + c(1, -1)*y.pad*diff(tmp.ylim);
## make final y coordinates by squeezing y.bin into tmp.ylim
if (is.null(grl.segs$y))
grl.segs$y = NA
if (all(is.na(grl.segs$y)))
{
grl.segs$y = affine.map(grl.segs$y.bin, tmp.ylim)
}
}
if (verbose) {
print('After agg')
print(Sys.time() - now)
} # xlim = c(min(window.segs$start), max(window.segs$end));
}
else
empty.plot = TRUE
#window.segs$end <- window.segs$end + 1 ##debug
## now collapse everything to 0, 1 based on windows.segs
winlim = range(c(window.segs$start, window.segs$end))
## xlim here is arbitrary, just needs to be > 1, helps us scale overlayed plots to window.segs boundaries (if many plots with different
## windows are drawn on the same canvas
grl.segs$pos1 = affine.map(grl.segs$pos1, winlim, ylim = xlim)
grl.segs$pos2 = affine.map(grl.segs$pos2, winlim, ylim = xlim)
window.segs$start = affine.map(window.segs$start, winlim, ylim = xlim)
window.segs$end = affine.map(window.segs$end, winlim, ylim = xlim)
}
else
empty.plot = TRUE
if (empty.plot)
{
if (is.null(windows))
stop('Either nonempty range data or windows must be provided')
mapped = gr.flatmap(GRanges(), windows, win.gap)
window.segs = mapped$window.segs
winlim = range(c(window.segs$start, window.segs$end))
window.segs$start = affine.map(window.segs$start, winlim, ylim = xlim)
window.segs$end = affine.map(window.segs$end, winlim, ylim = xlim)
# xlim = c(min(window.segs$start), max(window.segs$end));
if (is.null(ylim))
ylim = c(0, 1)
if (is.list(y) & is.null(ylim.subplot))
if (all(c('start', 'end') %in% names(y)))
ylim.subplot = c(y$start[1], y$end[1])
else
ylim.subplot = c(y[[1]], y[[2]])
if (is.null(ylim.subplot))
ylim.subplot = ylim
}
# if new plot will add (optional) axes and vertical lines separating windows
if (new.plot)
{
if (verbose) {
print('Before axis draw')
print(Sys.time() - now)
}
plot.blank(xlim = xlim, ylim = ylim, bg.col = bg.col);
new.axis = TRUE
}
if (is.na(sep.draw))
sep.draw = FALSE
if (sep.draw && length(windows)>1)
{
## rect(window.segs$end[1:(nrow(window.segs)-1)], rep(ylim[1], nrow(window.segs)-1),
## window.segs$start[2:(nrow(window.segs))], rep(ylim[2], nrow(window.segs)-1), border = 'white', col = sep.col)
if (any(width(windows)<=0))
warning('Some windows are width 0')
sep.loc = c(window.segs$start, window.segs$end);
if (!is.null(y.grid) && !all(is.na(y.grid)))
yrange = range(y.grid)
else if (!is.null(ylim.subplot))
yrange = ylim.subplot
else
yrange = range(grl.segs$y, na.rm = TRUE)
if (is.null(window.segs$border))
window.segs$border = 'white'
sep.x0 = window.segs$start[1:(nrow(window.segs))]
sep.x1 = window.segs$end[1:(nrow(window.segs))]
sep.y0 = rep(yrange[1], nrow(window.segs))
bgcol.l <- as.character(window.segs$col) ## JEREMIAH
bgborder.l <- as.character(window.segs$border) ## JEREMIAH
#rep(min(xaxis.pos.label, xaxis.pos, yrange[1]), nrow(window.segs))
sep.y1 = rep(yrange[2], nrow(window.segs))
rect(sep.x0, sep.y0, sep.x1, sep.y1, border = bgborder.l, col = bgcol.l) ## JEREMIAH added bgcol.l
segments(sep.x0, sep.y0, sep.x0, sep.y1, lty = sep.lty, lwd = sep.lwd)
segments(sep.x1, sep.y0, sep.x1, sep.y1, lty = sep.lty, lwd = sep.lwd)
## else
## {
## # rect(window.segs$start[1:(nrow(window.segs))], rep(ylim[1], nrow(window.segs)),
## #window.segs$end[1:(nrow(window.segs))], rep(ylim[2], nrow(window.segs)), border = 'white', col = sep.bg.col) ## MARCIN
## # window.segs$end[1:(nrow(window.segs))], rep(ylim[2], nrow(window.segs)), border = 'white', col = as.character(window.segs$col)) ## JEREMIAH
## # abline(v = sep.loc, col = 'gray', lty = sep.lty, lwd = sep.lwd);
## }
}
if (new.plot || new.axis)
{
if (is.null(xaxis.pos)) {
if (!is.null(ylim.subplot))
xaxis.pos = ylim.subplot[1]-0.05*diff(ylim.subplot)
else
xaxis.pos = ylim[1]+0.12*diff(ylim)
}
if (is.null(window.segs$col))
window.segs$col = sep.bg.col
if (is.null(xaxis.pos.label)) {
if (!is.null(ylim.subplot))
xaxis.pos.label = xaxis.pos - 0.04*diff(ylim.subplot)
else
xaxis.pos.label = xaxis.pos - 0.04*diff(ylim)
}
if (new.axis)
{
nwin = length(windows);
## draw the actual x axis
segments(window.segs$start, rep(xaxis.pos[1], nwin), window.segs$end, rep(xaxis.pos[1], nwin));
# if (!is.null(xaxis.suffix))
# if (is.na(xaxis.suffix) | nchar(xaxis.suffix)==0)
# xaxis.suffix = NULL
draw_x_ticks(xaxis.interval, windows, mapped, winlim, xlim, ylim, xaxis.pos, xaxis.suffix, xaxis.unit, xaxis.cex.tick, xaxis.ticklen, xaxis.round)
# then (label) text
newline <- ifelse(xaxis.newline, '\n', '')
width.text = ''
if (xaxis.width)
{
if (!is.null(xaxis.suffix))
width.text = paste('(', paste(prettyNum(ifelse(rep(xaxis.unit == 1, length(windows)),
width(windows), round(width(windows)/xaxis.unit, 2)), big.mark = ','), xaxis.suffix), ')', sep = '')
else
width.text = paste('(', prettyNum(ifelse(rep(xaxis.unit == 1, length(windows)),
width(windows), round(width(windows)/xaxis.unit, 2)), big.mark = ','), ')', sep = '')
}
begin.text = prettyNum(pmax(floor(1/xaxis.unit),
ifelse(rep(xaxis.unit == 1, length(windows)), start(windows), round(start(windows)/xaxis.unit, xaxis.round))),
big.mark = ',')
end.text = prettyNum(ifelse(rep(xaxis.unit == 1, length(windows)), end(windows),
round(end(windows)/xaxis.unit, xaxis.round)), big.mark = ',')
if (!xaxis.chronly) {
text(rowMeans(window.segs[, c('start', 'end')]), rep(xaxis.pos.label, nwin),
paste(xaxis.prefix, ' ', seqnames(windows), ':',newline,
begin.text,'-', newline,
end.text, ' ', xaxis.suffix, newline, width.text, sep = ''),
cex = xaxis.cex.label*0.7, srt = 0, adj = c(0.5, 0), srt=xaxis.label.angle)
} else {
text(rowMeans(window.segs[, c('start', 'end')]), rep(xaxis.pos.label, nwin),
paste(xaxis.prefix, ' ', seqnames(windows),
sep = ''),
cex = xaxis.cex.label*0.7, srt = 0, adj = c(0.5, 0), srt=xaxis.label.angle)
}
}
}
if (empty.plot)
{
if (verbose) {
print('Returning ..')
print(Sys.time() - now)
}
return(window.segs)
}
line.loc = NULL
if (!is.na(y.grid[1]))
{
if (is.logical(y.grid))
line.loc = pretty(ylim.grid,y.grid.pretty)
else if (length(y.grid)==1) ## only interval is specified
line.loc = seq(floor(ylim.grid[1]/y.grid)*y.grid, ceiling(ylim.grid[2]/y.grid)*y.grid, y.grid)
else ## specific grid lines are specified
line.loc = y.grid
if (is.null(names(line.loc)))
names(line.loc) = line.loc;
# abline(h = line.loc, col = y.grid.col, lty = y.grid.lty, lwd = y.grid.lwd)
segments(xlim[1], line.loc, xlim[2], line.loc, col = y.grid.col, lty = y.grid.lty, lwd = y.grid.lwd)
if (is.null(y.grid.cex))
y.grid.cex = NA
if (is.na(y.grid.cex))
y.grid.cex = 1
axis(2, at = line.loc, labels = names(line.loc), tick = TRUE, pos = line.loc[1], las = 2, cex.axis = y.grid.cex)
if (!is.null(ylab))
mtext(ylab, side = 2, at = mean(c(line.loc[1], line.loc[length(line.loc)])), line = 2, cex.lab = xaxis.cex.label)
}
if (length(grl)>0)
{
if (is.null(grl.segs$ywid))
grl.segs$ywid = 1
if (any(nix <- is.na(grl.segs$ywid)))
grl.segs$ywid[nix] = 1
if (!is.null(ylim.subplot))
tmp.ydiff = diff(ylim.subplot)*(1-2*y.pad)
else if (!is.null(line.loc))
tmp.ydiff = diff(range(line.loc))
else
tmp.ydiff = diff(range(grl.segs$y, na.rm = T))
if (tmp.ydiff==0)
tmp.ydiff = ylim
fact = 1.5*(1+length(unique(grl.segs$y)))
if (!is.null(line.loc))
fact = pmax(10, pmin(1000, fact))
## if ywid not provided then create from tmp.ydiff and then scale to new axis
if (is.null(ywid) || is.na(ywid))
{
ywid = tmp.ydiff / fact
if (!is.null(line.loc))
ywid = pmin(ywid, min(c(1, diff(sort(unique(grl.segs$y))))*2)) ## want to be able to see a minimal difference between data points
if (!is.null(line.loc)) ## don't want segments to be fatter than a grid unit
ywid = pmin(min(diff(line.loc))/2, ywid)
}
grl.segs$ywid = ywid * grl.segs$ywid
#########################################
## border / color logic finally
#########################################
if (is.null(gr.colorfield))
gr.colorfield = NA
if (gr.colorfield %in% names(grl.segs))
{
if (is.null(gr.colormap))
{
uval = unique(as.character(grl.segs[, gr.colorfield]))
gr.colormap = structure(alpha(brewer.master(length(uval)), 0.5), names = uval)
}
cols = gr.colormap[as.character(grl.segs[, gr.colorfield])];
grl.segs$col[!is.na(cols)] = cols[!is.na(cols)]
}
else
{
if (any(ix <- is.na(grl.segs$col)))
grl.segs$col[ix] = alpha('black', 0.5)
}
## border if unspecified will be a darker and less transparent version of the color
if (any(ix <- is.na(grl.segs$border)))
{
rgb = col2rgb(grl.segs$col[ix], alpha = TRUE)
grl.segs$border[ix] = rgb(rgb['red', ]/255, rgb['green', ]/255, rgb['blue', ]/255, alpha = 0.9)
}
if (leg.params$plot && length(gr.colormap)>0)
{
leg.params$x = leg.params$xpos * diff(xlim) + xlim[1]
leg.params$y = leg.params$ypos * diff(par('usr')[3:4]) + par('usr')[3]
leg.params$legend = names(gr.colormap)
if (circles) {
leg.params$col = gr.colormap
leg.params$pch = 16
}
else
leg.params$fill = gr.colormap
leg.params$border = gr.colormap
leg.params$xpos = leg.params$ypos = NULL
# if (length(gr.colormap)>legend.maxitems & legend.maxitems > 0)
# gr.colormap = gr.colormap[intersect(names(sort(table(grl.segs[, gr.colorfield]), decreasing = T)[1:legend.maxitems]),
# names(gr.colormap))]
do.call(graphics::legend, leg.params)
#if (circles)
# graphics::legend(legend.pos[1]*diff(xlim) + xlim[1], legend.pos[2]*diff(par('usr')[3:4]) + par('usr')[3], legend = names(gr.colormap), col = gr.colormap, border = gr.colormap, cex = legend.cex * 0.5, ncol = legend.ncol, xjust = legend.xjust, pch = 16, yjust = legend.yjust)
#else
# graphics::legend(legend.pos[1]*diff(xlim) + xlim[1], legend.pos[2]*diff(par('usr')[3:4]) + par('usr')[3], legend = names(gr.colormap), fill = gr.colormap, border = gr.colormap, cex = legend.cex * 0.5, ncol = legend.ncol, xjust = legend.xjust, yjust = legend.yjust)
}
if (draw.paths)
draw.backbone = FALSE
if (labels.suppress.gr)
grl.segs$label = NULL
if (!is.na(lines))
if (lines)
grl.segs = grl.segs[order(grl.segs$pos1), ]
draw.ranges(grl.segs, y = grl.segs$y, group = grl.segs$group, col = grl.segs$col, border = grl.segs$border, ylim = ylim, xlim = xlim, lwd = grl.segs$ywid, draw.backbone = draw.backbone, angle = angle, col.backbone = col.backbone, points = points, circles = circles, bars = bars, y0.bar = y0.bar, lines = lines, cex.label = gr.cex.label, srt.label = gr.srt.label, adj.label = gr.adj.label, ...)
## if draw.paths, will now draw connectors
##
if (draw.paths)
{
grl.segs = grl.segs[
order(grl.segs$grl.ix, # walk id
grl.segs$grl.iix, # walk segment
ifelse(grl.segs$strand=='+', 1, -1)*grl.segs$start) ## break ties in strand aware manner if walk segment spread across multiple (window) segments
, ]
ix.l = split(1:nrow(grl.segs), grl.segs$group)
grl.segs$ctype = NA; ## connector type
# keep track to see if next index is missing --> will allow us to draw dotted connectors for these ..
# "missing" indices occur when we are focusing on windows and potentially removing some of the
# ranges in the contig
grl.segs$next.missing = !((grl.segs$query.id+1) %in% grl.segs$query.id) & !grl.segs$last
grl.segs$prev.missing = !((grl.segs$query.id-1) %in% grl.segs$query.id) & !grl.segs$first
if (is.null(grl.segs$is.cycle))
grl.segs$is.cycle = FALSE;
connector.args = do.call('rbind', lapply(ix.l, function(ix)
{
# find runs where start[i+1]>end[i] and strand[i] == strand[i+1] = '+'
# and end[i+1]<start[i] and strand[i] == strand[i+1] = '-'
out = NULL;
if (length(ix)>1)
{
out = data.frame(ix0 = ix[-length(ix)], ix1 = ix[-1], type = 'U', sign = '+', cyclic = F, stringsAsFactors = F);
discordant = grl.segs$y.bin[ix[-length(ix)]] != grl.segs$y.bin[ix[-1]]
out$type[discordant & grl.segs$strand[ix[-1]] == grl.segs$strand[ix[-length(ix)]]] = 'S'
out$type[discordant & grl.segs$strand[ix[-1]] != grl.segs$strand[ix[-length(ix)]]] = 'S'
}
if (grl.segs$next.missing[ix[length(ix)]]) ## make bridge to nowhere
out = rbind(out, data.frame(ix0 = ix[length(ix)], ix1 = NA, type = 'S', cyclic = F,
sign = grl.segs$strand[ix[length(ix)]], stringsAsFactors = F))
if (grl.segs$prev.missing[ix[1]]) ## make bridge from nowhere
out = rbind(out, data.frame(ix0 = NA, ix1 = ix[1], type = 'S', cyclic = F,
sign = grl.segs$strand[ix[1]], stringsAsFactors = F))
if (grl.segs$is.cycle[ix[length(ix)]])
{
if (grl.segs$y.bin[ix[length(ix)]] == grl.segs$y.bin[ix[1]])
out = rbind(out, data.frame(ix0 = ix[length(ix)], ix1 = ix[1], type = 'U', cyclic = T, sign = '-', stringsAsFactors = F))
else if (grl.segs$strand[ix[length(ix)]] == '-' & grl.segs$strand[ix[1]] == '+')
out = rbind(out, data.frame(ix0 = ix[length(ix)], ix1 = ix[1], type = 'S', cyclic = T, sign = '+', stringsAsFactors = F))
else if (grl.segs$strand[ix[length(ix)]] == '+' & grl.segs$strand[ix[1]] == '-')
out = rbind(out, data.frame(ix0 = ix[length(ix)], ix1 = ix[1], type = 'S', cyclic = T, sign = '-', stringsAsFactors = F))
else if (grl.segs$strand[ix[length(ix)]] == '-' & grl.segs$strand[ix[1]] == '-')
out = rbind(out, data.frame(ix0 = ix[length(ix)], ix1 = ix[1], type = 'S', cyclic = T, sign = '+', stringsAsFactors = F))
else if (grl.segs$strand[ix[length(ix)]] == '+' & grl.segs$strand[ix[1]] == '+')
out = rbind(out, data.frame(ix0 = ix[length(ix)], ix1 = ix[1], type = 'S', cyclic = T, sign = '-', stringsAsFactors = F))
}
return(out)
}))
if (!is.null(connector.args))
{
path.h = path.cex.h * rep(diff(par('usr')[1:2])/50, nrow(connector.args))
if (any(connector.args$type == 'S'))
path.h[connector.args$type == 'S'] = 2*path.h[connector.args$type == 'S']
path.v = rep(path.cex.v, nrow(connector.args))
path.v[is.na(connector.args$ix0) | is.na(connector.args$ix1)] = path.cex.v*2*grl.segs$ywid[connector.args$ix0[is.na(connector.args$ix0) | is.na(connector.args$ix1)]]
# path.v[is.na(connector.args$ix0) | is.na(connector.args$ix1)] = path.cex.v*2*grl.segs$ywid[connector.args$ix0[is.na(connector.args$ix0) | is.na(connector.args$ix1)]]
# path.v[!is.na(connector.args$ix0)] = path.cex.v*2*grl.segs$ywid[connector.args$ix0[!is.na(connector.args$ix0)]]
# path.v[!is.na(connector.args$ix1)] = pmax(path.v[!is.na(connector.args$ix1)],
# path.cex.v*2*grl.segs$ywid[connector.args$ix1[!is.na(connector.args$ix1)]])
connector.args$path.col = path.col
connector.args$path.lwd = path.lwd
## override path.col and path.lwd if provided in grl.segs
if (!is.null(grl.segs$path.col))
connector.args$path.col[!is.na(connector.args$ix0)] = grl.segs$path.col[connector.args$ix0[!is.na(connector.args$ix0)]]
if (!is.null(grl.segs$path.lwd))
connector.args$path.lwd[!is.na(connector.args$ix0)] = grl.segs$path.lwd[connector.args$ix0[!is.na(connector.args$ix0)]]
connector.args$strand0[!is.na(connector.args$ix0)] = grl.segs$strand[connector.args$ix0[!is.na(connector.args$ix0)]]
connector.args$strand1[!is.na(connector.args$ix1)] = grl.segs$strand[connector.args$ix1[!is.na(connector.args$ix1)]]
connector.args$strand0[is.na(connector.args$ix0)] = grl.segs$strand[connector.args$ix1[is.na(connector.args$ix0)]]
connector.args$strand1[is.na(connector.args$ix1)] = grl.segs$strand[connector.args$ix0[is.na(connector.args$ix1)]]
connector.args$strand0 = ifelse(connector.args$strand0 == '*', '+', connector.args$strand0)
connector.args$strand1 = ifelse(connector.args$strand1 == '*', '+', connector.args$strand1)
connector.args$x0[!is.na(connector.args$ix0) & connector.args$strand0=='+'] =
grl.segs$pos2[connector.args$ix0[!is.na(connector.args$ix0) & connector.args$strand0=='+']]
connector.args$x0[!is.na(connector.args$ix0) & connector.args$strand0=='-'] =
grl.segs$pos1[connector.args$ix0[!is.na(connector.args$ix0) & connector.args$strand0=='-']]
connector.args$x1[!is.na(connector.args$ix1) & connector.args$strand1=='+'] =
grl.segs$pos1[connector.args$ix1[!is.na(connector.args$ix1) & connector.args$strand1=='+']]
connector.args$x1[!is.na(connector.args$ix1) & connector.args$strand1=='-'] =
grl.segs$pos2[connector.args$ix1[!is.na(connector.args$ix1) & connector.args$strand1=='-']]
# bridges from nowhere
connector.args$x0[is.na(connector.args$ix0) & connector.args$strand1=='+'] =
grl.segs$pos1[connector.args$ix1[is.na(connector.args$ix0) & connector.args$strand1=='+']] - path.h[is.na(connector.args$ix0) & connector.args$strand1=='+']
connector.args$x0[is.na(connector.args$ix0) & connector.args$strand1=='-'] =
grl.segs$pos2[connector.args$ix1[is.na(connector.args$ix0) & connector.args$strand1=='-']] + path.h[is.na(connector.args$ix0) & connector.args$strand1=='-']
# bridges to nowhere
connector.args$x1[is.na(connector.args$ix1) & connector.args$strand0=='+'] =
grl.segs$pos2[connector.args$ix0[is.na(connector.args$ix1) & connector.args$strand0=='+']] + path.h[is.na(connector.args$ix1) & connector.args$strand0=='+']
connector.args$x1[is.na(connector.args$ix1) & connector.args$strand0=='-'] =
grl.segs$pos1[connector.args$ix0[is.na(connector.args$ix1) & connector.args$strand0=='-']] - path.h[is.na(connector.args$ix1) & connector.args$strand0=='-']
connector.args$y0[!is.na(connector.args$ix0)] = grl.segs$y[connector.args$ix0[!is.na(connector.args$ix0)]]
connector.args$y1[!is.na(connector.args$ix1)] = grl.segs$y[connector.args$ix1[!is.na(connector.args$ix1)]]
connector.args$y0[is.na(connector.args$ix0)] = grl.segs$y[connector.args$ix1[is.na(connector.args$ix0)]] - 0.25*path.v[is.na(connector.args$ix0)]
connector.args$y1[is.na(connector.args$ix1)] = grl.segs$y[connector.args$ix0[is.na(connector.args$ix1)]] + 0.25*path.v[is.na(connector.args$ix1)]
connector.args$lty = 1;
connector.args$lty[grl.segs$next.missing[connector.args$ix0[!is.na(connector.args$ix0)]] &
!connector.args$cyclic[!is.na(connector.args$ix0)]] = 3
connector.args$lty[grl.segs$prev.missing[connector.args$ix1[!is.na(connector.args$ix1)]] &
connector.args$cyclic[!is.na(connector.args$ix1)]] = 3
connector.args$lty[is.na(connector.args$ix0) | is.na(connector.args$ix1)] = 3; ## label all bridge to / from nowhere with dotted line
##
path.v[connector.args$y0 == connector.args$y1] = 0
path.h[connector.args$y0 == connector.args$y1] = 0
path.v[connector.args$cyclic] = path.cex.v*2*grl.segs$ywid[connector.args$ix1[connector.args$cyclic]]
path.h[connector.args$cyclic] = path.cex.h * diff(par('usr')[1:2])/20
## workaround current lines() limitation in connectors
if (any(lty3 <- connector.args$lty == 3))
connectors(connector.args$x0[lty3], connector.args$y0[lty3], connector.args$strand0[lty3],
connector.args$x1[lty3], connector.args$y1[lty3], ifelse(connector.args$strand1[lty3] == '+', '-', '+'),
type = connector.args$type[lty3],
h = path.h[lty3], v = path.v[lty3],
lwd = connector.args$path.lwd[lty3],
lty = connector.args$lty[lty3],
col = connector.args$path.col[lty3],
col.arrow = path.col.arrow,
cex.arrow = grl.segs$ywid[1]*path.cex.arrow, f.arrow = T)
if (any(lty1 <- connector.args$lty == 1))
connectors(connector.args$x0[lty1], connector.args$y0[lty1], connector.args$strand0[lty1],
connector.args$x1[lty1], connector.args$y1[lty1], ifelse(connector.args$strand1[lty1] == '+', '-', '+'),
type = connector.args$type[lty1],
h = path.h[lty1], v = path.v[lty1],
lty = connector.args$lty[lty1],
lwd = connector.args$path.lwd[lty1],
col = connector.args$path.col[lty1],
col.arrow = path.col.arrow,
cex.arrow = grl.segs$ywid[1]*path.cex.arrow, f.arrow = T)
#text(connector.args$x1, connector.args$y1, paste(connector.args$ix0, ' (', grl.segs$group[connector.args$ix0], '), ', connector.args$ix1, ' (', grl.segs$group[connector.args$ix1], '), ', connector.args$type, ' ', connector.args$sign, sep = ''))
# text(connector.args$x1, connector.args$y1-path.v, paste(grl.segs$group[connector.args$ix0], connector.args$type, ' ', connector.args$sign, sep = ''), )
}
}
if (!is.null(edges))
if (nrow(edges)>0 & all(c('from', 'to') %in% colnames(edges)))
{
if (data.table::is.data.table(edges))
edges = as.data.frame(edges)
if (is.null(edges$col))
edges$col = NA
if (is.null(edges$lty))
edges$lty = NA
if (is.null(edges$lwd))
edges$lwd = NA
if (is.null(edges$h))
edges$h = 1
if (is.null(edges$cex.arrow))
edges$cex.arrow = 1
edges.og = edges
edges = edges[edges$to %in% grl.segs$group | edges$from %in% grl.segs$group, ]
if (nrow(edges)>0)
{
if (any(ix <- is.na(edges$col)))
edges$col[ix] = 'black'
if (any(ix <- is.na(edges$lwd)))
edges$lwd[ix] = 1
if (any(ix <- is.na(edges$lty)))
edges$lty[ix] = 1
## now translate $from $to from group to gr id, and choosing last range in group for "from" indices
## and first range in group for "to" indices putting NA's when gr not in grl.segs
first.ix = which(grl.segs$first)
last.ix = which(grl.segs$last)
grl.segs$to.gr = grl.segs$from.gr = 1:nrow(grl.segs)
edges$edge.id = 1:nrow(edges)
edges = merge(merge(edges, grl.segs[last.ix, c('group', 'from.gr')], by.x = 'from', by.y = 'group', all.x = T),
grl.segs[first.ix, c('group', 'to.gr')], by.x = 'to', by.y = 'group', all.x = T)
# edges$to.gr = first.ix[match(edges$to, grl.segs[first.ix, ]$group)]
# edges$from.gr = last.ix[match(edges$from, grl.segs[last.ix, ]$group)]
## replicate edges that connect input gr's that have multiple instantiations
## figure out which grs are clipped
grl.segs$clipped.start = !(start(gr)[grl.segs$query.id] == grl.segs$start)
grl.segs$clipped.end = !(end(gr)[grl.segs$query.id] == grl.segs$end)
clipped.to = ((grl.segs$clipped.start[edges$to.gr] & grl.segs$strand[edges$to.gr] == '+') |
(grl.segs$clipped.end[edges$to.gr] & grl.segs$strand[edges$to.gr] == '-'))
clipped.from = ((grl.segs$clipped.start[edges$from.gr] & grl.segs$strand[edges$from.gr] == '-') |
(grl.segs$clipped.end[edges$from.gr] & grl.segs$strand[edges$from.gr] == '+'))
clipped.to[is.na(clipped.to)] = F
clipped.from[is.na(clipped.from)] = F
## now determine start and end points based on signs of connectors
to.ix = !is.na(edges$to.gr) & !clipped.to
from.ix = !is.na(edges$from.gr) & !clipped.from
edges$x.pos.from = edges$x.pos.to = edges$y.pos.from = edges$y.pos.to = edges$dir.from = edges$dir.to = NA
# from.ix = !is.na(edges$fromx.gr)
# to.ix = !is.na(edges$to.gr)
if (any(to.ix)) ## connect to left end if + and right end if -
{
edges$x.pos.to[to.ix] = ifelse(grl.segs$strand[edges$to.gr[to.ix]] != '-',
grl.segs$pos1[edges$to.gr[to.ix]], grl.segs$pos2[edges$to.gr[to.ix]])
edges$y.pos.to[to.ix] = grl.segs$y[edges$to.gr[to.ix]]
edges$dir.to[to.ix] = ifelse(grl.segs$strand[edges$to.gr[to.ix]] != '-',
'-', '+')
}
if (any(from.ix)) ## connect from right end if + and left end if -
{
edges$x.pos.from[from.ix] = ifelse(grl.segs$strand[edges$from.gr[from.ix]] != '-',
grl.segs$pos2[edges$from.gr[from.ix]], grl.segs$pos1[edges$from.gr[from.ix]])
edges$y.pos.from[from.ix] = grl.segs$y[edges$from.gr[from.ix]]
edges$dir.from[from.ix] = ifelse(grl.segs$strand[edges$from.gr[from.ix]] != '-',
'+', '-')
}
## now let NA endpoints "dangle"
dangle.w = diff(par('usr')[1:2])/100
edges$dangle = F
if (is.null(edges$dangle.w))
edges$dangle.w = 1
edges$dangle.w = edges$dangle.w * dangle.w
if (any(!from.ix))
{
edges$x.pos.from[!from.ix] =
edges$x.pos.to[!from.ix] + ifelse(grl.segs$strand[edges$to.gr[!from.ix]] != '-', -1, 1)*edges$dangle.w[!from.ix]
edges$y.pos.from[!from.ix] = edges$y.pos.to[!from.ix]
edges$from.gr[!from.ix] = edges$from.gr[!from.ix]
edges$dir.from[!from.ix] = ifelse(edges$dir.to[!from.ix] != '-', '-', '+')
edges$dangle[!from.ix] = T
}
if (any(!to.ix))
{
edges$x.pos.to[!to.ix] =
edges$x.pos.from[!to.ix] + ifelse(grl.segs$strand[edges$from.gr[!to.ix]] != '-', 1, -1)*edges$dangle.w[!to.ix]
edges$y.pos.to[!to.ix] = edges$y.pos.from[!to.ix]
edges$to.gr[!to.ix] = edges$to.gr[!to.ix]
edges$dir.to[!to.ix] = ifelse(edges$dir.from[!to.ix] != '-', '-', '+')
edges$dangle[!to.ix] = T
}
edges$from.ix = from.ix
edges$to.ix = to.ix
edges = edges[order(edges$dangle), ]; ## hack to make sure we prefer non-dangling versions of edges when we have a choice (TODO: MAKE NON HACKY)
dup = duplicated(cbind(edges$edge.id, edges$from.gr)) | duplicated(cbind(edges$edge.id, edges$to.gr))
edges = edges[(edges$from.ix | edges$to.ix) & !dup, ]
if (nrow(edges)>0)
{
if (!is.null(ylim.subplot))
tmp.ydiff = diff(ylim.subplot)*(1-2*y.pad)
else
tmp.ydiff = diff(range(grl.segs$y, na.rm = T))
# edges$type = ifelse(edges$y.pos.from == edges$y.pos.to, 'U', 'S')
uthresh = max(grl.segs$ywid)
edges$type = ifelse(abs(edges$y.pos.from - edges$y.pos.to) <= uthresh, 'U', 'S')
edges$f.arrow = T
edges$b.arrow = F
if (is.null(edges$v))
edges$v = 1
edges$v = ifelse(abs(edges$y.pos.from - edges$y.pos.to)<=uthresh,
2*uthresh*affine.map(abs(edges$x.pos.to-edges$x.pos.from), xlim = c(0, diff(par('usr')[1:2])), ylim = c(0.5, 1.0)),
abs(edges$y.pos.from-edges$y.pos.to)/2)*edges$v
# edges$v[is.na(edges$v)] = 0
edges$h = dangle.w*edges$h * affine.map(abs(edges$x.pos.to-edges$x.pos.from), xlim = c(0, diff(par('usr')[1:2])), ylim = c(0.5, 1.0))
make.flat = edges$type == 'U' & edges$dir.to != edges$dir.from
if (!is.null(edges$not.flat)) ## allow edges specified as $not.flat to have height, unless dangling
if (length(which(ix <- edges$not.flat & !edges$dangle))>0)
make.flat[ix] = F
if (any( ix <- make.flat ))
{
edges$v[ix] = 0
edges$h[ix] = 0
}
if (!is.null(edges$v.override))
if (any(ix <- (!is.na(edges$v.override) & edges$to.ix & edges$from.ix)))
edges$v[ix] = edges$v.override[ix]
if (!is.null(edges$h.override))
if (any(ix <- (!is.na(edges$h.override) & edges$to.ix & edges$from.ix)))
edges$h[ix] = dangle.w*edges$h.override[ix]
connectors(edges$x.pos.from, edges$y.pos.from, edges$dir.from,
edges$x.pos.to, edges$y.pos.to, edges$dir.to,
v = edges$v, h = edges$h, type = edges$type,
f.arrow = T, b.arrow = F,
cex.arrow = 0.2*edges$cex.arrow,
col.arrow = edges$col,
lwd = edges$lwd, lty = edges$lty, col = edges$col)
if (!is.null(edges$label))
{
text((edges$x.pos.from + edges$x.pos.to)/2, edges$y.pos.from + 0.5*edges$v*sign(edges$y.pos.to - edges$y.pos.from + 0.01),
edges$label, adj = c(0.5, 0.5), col = 'black')
}
}
}
}
if (verbose) {
print('After draw')
print(Sys.time() - now)
}
if (nrow(grl.segs)>0 & !is.null(grl.props$grl.labels) & !labels.suppress.grl)
{
if (!draw.paths)
{
pos1 = vaggregate(pos1 ~ group, data = grl.segs, FUN = min);
pos2 = vaggregate(pos2 ~ group, data = grl.segs, FUN = max);
ywid = vaggregate(ywid ~ group, data = grl.segs, FUN = max);
y = vaggregate(y ~ group, data = grl.segs, FUN = mean);
grl.segs.u = data.frame(group = names(pos1), pos1, pos2, y, ywid);
grl.segs.u$grl.labels = grl.props$grl.labels[match(grl.segs.u$group, grl.props$group)]
grl.segs.u$grl.cols = grl.props$grl.cols[match(grl.segs.u$group, grl.props$group)]
grl.segs.u$grl.cex = grl.props$grl.cex[match(grl.segs.u$group, grl.props$group)]
rownames(grl.segs.u) = grl.segs.u$group;
if (adj.label[2] == 0)
text(grl.segs.u$pos1, grl.segs.u$y+grl.segs.u$ywid + 0.005*diff(ylim),
grl.segs.u$grl.labels, adj = adj.label, cex = grl.segs.u$grl.cex, col = grl.segs.u$grl.cols)
if (adj.label[2] == 0.5)
text(grl.segs.u$pos1-0.005*diff(xlim), grl.segs.u$y,
grl.segs.u$grl.labels, adj = adj.label, cex = grl.segs.u$grl.cex, col = grl.segs.u$grl.cols)
else
text(grl.segs.u$pos1, grl.segs.u$y-grl.segs.u$ywid - 0.005*diff(ylim),
grl.segs.u$grl.labels, adj = adj.label, cex = grl.segs.u$grl.cex, col = grl.segs.u$grl.cols)
}
else
{
pos1 = vaggregate(pos1 ~ group, data = grl.segs, FUN = min);
pos2 = vaggregate(pos2 ~ group, data = grl.segs, FUN = max);
y0 = vaggregate(y ~ group, data = grl.segs, FUN = min);
y1 = vaggregate(y ~ group, data = grl.segs, FUN = max);
grl.segs.u = data.frame(group = names(pos1), pos1, pos2, y0, y1);
grl.segs.u$grl.labels = grl.props$grl.labels[match(grl.segs.u$group, grl.props$group)]
grl.segs.u$grl.cols = grl.props$grl.cols[match(grl.segs.u$group, grl.props$group)]
text(adj.label[1]*grl.segs.u$pos1 + (1-adj.label[1])*grl.segs.u$pos2, adj.label[2]*grl.segs.u$y0 + (1-adj.label[2])*grl.segs.u$y1,
grl.segs.u$grl.labels, adj = adj.label, cex = grl.segs.u$grl.cex, col = grl.segs.u$grl.cols)
}
}
}
return(window.segs)
}
#####################################
#' @name barbs
#' @title barbs
#' @description
#'
#' (improvement over pencils)
#' used to represent directional ranges. Shape consists of two horizontally reflected parallelograms
#' with a specified angle. The middle segment of the parallelogram represents the provided range.
#'
#' Parameterized by x0, x1 (range values), y0, y1 (thickness), angle, where angle is between -90 and 90
#'
#' eg angle of 0 = rectangle
#' angle of 45 = rightward pointing barb
#' angle of -45 = leftward pointing barb
#'
#' Note: the shape will NOT be necessarily contained inside the box outlined by x0, x1, y0, y1
#' @author Marcin Imielinski
#' @keywords internal
barbs = function(x0, y0, x1, y1, angle = 0, ...)
{
y.thick = (y1-y0)/2;
y.midpoint = rowMeans(cbind(y0, y1))
angle = pmax(-80, pmin(80, angle))/180*pi;
# to make the angle visible irrespective of plot coordinates
# we need to compute the x.offset in viewer coordinates (ie inches)
x.wid = diff(par('usr')[1:2])
y.wid = diff(par('usr')[3:4])
x.offset = tan(angle)*y.thick/y.wid*par('pin')[2]*x.wid/par('pin')[1]
x.vertices = rbind(x0, x0-x.offset, x1-x.offset, x1, x1-x.offset, x0-x.offset, x0, NA)
y.vertices = rbind(y.midpoint, y1, y1, y.midpoint, y0, y0, y.midpoint, NA);
if (ncol(x.vertices)==1)
x.vertices = rep(x.vertices, ncol(y.vertices))
else if (ncol(y.vertices)==1)
y.vertices = rep(y.vertices, ncol(x.vertices))
## convert lwd.border 0 to NA. JEREMIAH
other.args <- list(...)
if (!any(other.args$lwd > 0, na.rm = TRUE)) {
other.args$border = NA
other.args[['lwd']] <- NULL
}
if (!is.null(other.args$lwd))
{
if (any(ix <- is.na(other.args$lwd)))
other.args$lwd[ix] = 1
}
do.call('polygon', c(list(x=as.numeric(x.vertices), y=as.numeric(y.vertices)), other.args[names(other.args) %in% c("lwd", "border", "col")]))
}
#' @name bernsteinp
#' @title bernsteinp
#' @description
#' bernsteinp
#'
#' computes matrix of bernstein polynomials of degree n at m points t in the
#' interval 0, 1
#'
#' out[i,j] = choose(n, j)*t[i]^j*(1-t[i])^(n-j)
#' where t = seq(0, 1, length.out = m) and j in 0 .. n
#'
#' this matrix can be multiplied by a n x 2 matrix of control points to yield
#' an m x 2 matrix of m equally spaced points along the curve
#' @author Marcin Imielinski
#' @keywords internal
bernsteinp = function(n, m)
{
m = round(m);
if (m<0)
return(NULL)
t = seq(0, 1, length.out = m)
out = matrix(sapply(0:(n-1), function(i) choose(n-1,i)*t^i*(1-t)^(n-1-i)), nrow = m);
return(out)
}
#' @name connectors
#' @title connectors
#' @description
#'
#' connectors
#'
#' draws bezier connectors between pairs of signed points points a = (x0, y0, s0) and b(x1, y1, s1)
#'
#' "S" connectors will use bezier control points that are located in an intermediate y location
#' between y0 and y1, and "U" connectors will use bezier control points that are located in an
#' intermediate x location between x0 and x1.
#'
#' type can be "U" or "S", if "S", then sign of corresponding v value is ignored (ie control points will always
#' be vertically in between a and b, in addition, neither intermeidate control will never be more than 1/3 of the y distance
#' between a and b
#'
#' The sign of the endpoints refers to the direction from which the connector approaches the point, s0 = 1 (or '+') means from the right
#' and s0 = 0 means from the left, similarly for s1
#'
#' h = horizonal "bulge" of signed connector, ie depending on the sign the extent to which the connector will extend to the left
#' or the right of points an and b
#' v = vertical bulge (only applies to U connectors), how far above max(y0, y1) the U connector will bulge if v>0 or how far below min(y0, y1)
#' the connector will bulge.
#'
#' nsteps determines how many segments used to approximate the curve
#'
#' all args except nsteps can be vecorized
#'
#' @author Marcin Imielinski
#' @keywords internal
connectors = function(x0, y0, s0 = 1, x1, y1, s1 = 1, v = 0.1, h = 0.1,
type = "S", f.arrow = T, b.arrow = F, nsteps = 100,
cex.arrow = 1, col.arrow = 'black', lwd = 1, lty = 1, col = 'black')
{
B = rbind(bernsteinp(6, nsteps), NA); # bernstein basis for bezier with 6 control points
dat = data.frame(x0, y0, s0, x1, y1, s1, type, v, h, f.arrow, b.arrow, cex.arrow, lwd, lty, stringsAsFactors = F, col = col, col.arrow = col.arrow)
cpoints.names =
matrix(paste('c', rep(1:ncol(B), each = 2), c('.x', '.y'), sep = ''), nrow = 2, ncol = ncol(B)-2, dimnames = list(c('x', 'y'), 1:(ncol(B)-2)))
dat[, as.vector(cpoints.names)] = NA
if (length(dat) == 0)
return(NULL)
scon = dat$type == "S"
ucon = dat$type == "U"
if (is.numeric(dat$s0))
dat$s0 = sign(dat$s0)
else
dat$s0 = sign((dat$s0 == '+')-0.5)
if (is.numeric(dat$s1))
dat$s1 = sign(dat$s1)
else
dat$s1 = sign((dat$s1 == '+')-0.5)
if (any(scon))
{
dat.s = dat[scon, , drop = FALSE]
dat.s$v = pmin(abs(dat.s$v), abs(dat.s$y1-dat.s$y0)/2)
dat.s$h = abs(dat.s$h)
dat.s[, cpoints.names['x', 1]] = dat.s$x0 + dat.s$s0*dat.s$h
dat.s[, cpoints.names['y', 1]] = dat.s$y0
## middle control points will form an "S" or "C" depending on s0 and s1
dat.s[, cpoints.names['x', 2]] = dat.s[, cpoints.names['x', 1]]
dat.s[, cpoints.names['y', 2]] = dat.s[, cpoints.names['y', 1]] + ifelse(dat.s$y0>dat.s$y1, -1, 1)*dat.s$v
dat.s[, cpoints.names['x', 4]] = dat.s$x1 + dat.s$s1*dat.s$h
dat.s[, cpoints.names['y', 4]] = dat.s$y1
dat.s[, cpoints.names['x', 3]] = dat.s[, cpoints.names['x', 4]]
dat.s[, cpoints.names['y', 3]] = dat.s[, cpoints.names['y', 4]] + ifelse(dat.s$y0>dat.s$y1, 1, -1)*dat.s$v
dat[scon, ] = dat.s
}
if (any(ucon))
{
dat.u = dat[ucon, , drop = FALSE]
# dat.u$h = pmin(abs(dat.u$h), abs(dat.u$x1-dat.u$x0)/3)
dat.u$h = abs(dat.u$h)
dat.u[, cpoints.names['x', 1]] = ifelse(dat.u$x0 < dat.u$x1, pmin((dat.u$x0 + dat.u$x1)/2, dat.u$x0 + dat.u$s0*dat.u$h), pmax((dat.u$x0 + dat.u$x1)/2, dat.u$x0 + dat.u$s0*dat.u$h))
# dat.u[, cpoints.names['x', 1]] = dat.u$x0 + dat.u$s0*dat.u$h
dat.u[, cpoints.names['y', 1]] = dat.u$y0
## middle control points will form verticle bubble portion of U (if v<0) or upside down U (if v>0)
dat.u[, cpoints.names['x', 2]] = dat.u$x0
dat.u[, cpoints.names['y', 2]] = ifelse(dat.u$v>=0, pmax(dat.u$y0, dat.u$y1) + dat.u$v, pmin(dat.u$y0, dat.u$y1) + dat.u$v)
dat.u[, cpoints.names['x', 3]] = dat.u$x1
dat.u[, cpoints.names['y', 3]] = dat.u[, cpoints.names['y', 2]]
dat.u[, cpoints.names['x', 4]] = ifelse(dat.u$x0 < dat.u$x1, pmax((dat.u$x0 + dat.u$x1)/2, dat.u$x1 + dat.u$s1*dat.u$h), pmin((dat.u$x0 + dat.u$x1)/2, dat.u$x1 + dat.u$s1*dat.u$h))
# dat.u[, cpoints.names['x', 4]] = dat.u$x1 + dat.u$s1*dat.u$h
dat.u[, cpoints.names['y', 4]] = dat.u$y1
dat[ucon, ] = dat.u
}
## draw bezier around control points .. draw separately for each unique combo of lwd, lty, col
formats = factor(paste(dat$lwd, dat$lty, dat$col))
lapply(levels(formats), function(lev)
{
lev.ix = which(formats == lev)
x = as.vector(t(as.matrix(dat[lev.ix, c('x0', cpoints.names['x', ], 'x1')])))
y = as.vector(t(as.matrix(dat[lev.ix, c('y0', cpoints.names['y', ], 'y1')])))
ix = seq(1, length(x), ncol(B))
XY = do.call('rbind', lapply(ix, function(i) B %*% cbind(x[i:(i+ncol(B)-1)], y[i:(i+ncol(B)-1)])))
lines(XY, lwd = dat$lwd[lev.ix], lty = dat$lty[lev.ix], col = as.character(dat$col)[lev.ix])
})
x.wid = diff(par('usr')[1:2])
y.wid = diff(par('usr')[3:4])
arrow.height = cex.arrow*0.05*y.wid
x.offset = tan(120)*arrow.height/y.wid*par('pin')[2]*x.wid/par('pin')[1]
## add forward arrow heads if any
if (any(dat$f.arrow))
{
tmp = dat$x1[dat$f.arrow] + dat$s1[dat$f.arrow]*x.offset
xcoord = rbind(dat$x1[dat$f.arrow], tmp, tmp, NA);
ycoord = rbind(dat$y1[dat$f.arrow], dat$y1[dat$f.arrow] + arrow.height/2, dat$y1[dat$f.arrow] - arrow.height/2, NA);
polygon(xcoord, ycoord, col = as.character(dat$col.arrow[dat$f.arrow]), border = as.character(dat$col.arrow[dat$f.arrow]))
}
if (any(dat$b.arrow))
{
tmp = dat$x0[dat$b.arrow] + dat$s0[dat$b.arrow]*x.offset
xcoord = rbind(dat$x0[dat$b.arrow], tmp, tmp, NA);
ycoord = rbind(dat$y0[dat$b.arrow], dat$y0[dat$b.arrow] + arrow.height/2, dat$y0[dat$b.arrow] - arrow.height/2, NA);
polygon(xcoord, ycoord, col = as.character(dat$col.arrow[dat$b.arrow]), border = as.character(dat$col.arrow[dat$b.arrow]))
}
}
#' @name affine.map
#' @title affine.map
#' @description
#'
#'
#' affinely maps 1D points in vector x from interval xlim to interval ylim,
#' ie takes points that lie in
#' interval xlim and mapping onto interval ylim using linear / affine map defined by:
#' (x0,y0) = c(xlim(1), ylim(1)),
#' (x1,y1) = c(xlim(2), ylim(2))
#' (using two point formula for line)
#' useful for plotting.
#'
#' if cap.max or cap.min == T then values outside of the range will be capped at min or max
#' @rdname affine-map-methods
#' @author Marcin Imielinski
#' @keywords internal
affine.map = function(x, ylim = c(0,1), xlim = c(min(x), max(x)), cap = F, cap.min = cap, cap.max = cap, clip = T, clip.min = clip, clip.max = clip)
{
# xlim[2] = max(xlim);
# ylim[2] = max(ylim);
if (xlim[2]==xlim[1])
y = rep(mean(ylim), length(x))
else
y = (ylim[2]-ylim[1]) / (xlim[2]-xlim[1])*(x-xlim[1]) + ylim[1]
if (cap.min)
y[x<min(xlim)] = ylim[which.min(xlim)]
else if (clip.min)
y[x<min(xlim)] = NA;
if (cap.max)
y[x>max(xlim)] = ylim[which.max(xlim)]
else if (clip.max)
y[x>max(xlim)] = NA;
return(y)
}
#' @name alpha
#' @title alpha
#' @description
#' Give transparency value to colors
#'
#' Takes provided colors and gives them the specified alpha (ie transparency) value
#'
#' @author Marcin Imielinski
#' @param col RGB color
#' @keywords internal
alpha = function(col, alpha)
{
col.rgb = col2rgb(col)
out = rgb(red = col.rgb['red', ]/255, green = col.rgb['green', ]/255, blue = col.rgb['blue', ]/255, alpha = alpha)
names(out) = names(col)
return(out)
}
#' Blends colors
#'
#' @param cols colors to blend
#' @keywords internal
blend = function(cols)
{
col.rgb = rowMeans(col2rgb(cols))
out = rgb(red = col.rgb['red']/255, green = col.rgb['green']/255, blue = col.rgb['blue']/255)
return(out)
}
#' @name col.scale
#' @title col.scale
#' @description
#'
#' Assign rgb colors to numeric data
#'
#' Assigns rgb colors to numeric data values in vector "x".. maps scalar values
#' in val.range (default c(0,1)) to a linear color scale of between col.min (default white)
#' and col.max (default black), each which are length 3 vectors or characters. RGB values are scaled between 0 and 1.
#' Values below and above val.min and val.max are mapped to col.max and col.max respectively
#'
#' @author Marcin Imielinski
#' @keywords internal
col.scale = function(x, val.range = c(0, 1), col.min = 'white', col.max = 'black', na.col = 'white',
invert = F # if T flips rgb.min and rgb.max
)
{
## NOTE fix
error = NULL
if (!is.numeric(col.min))
if (is.character(col.min))
col.min = col2rgb(col.min)/255
else
error('Color should be either length 3 vector or character')
if (!is.numeric(col.max))
if (is.character(col.max))
col.max = col2rgb(col.max)/255
else
error('Color should be either length 3 vector or character')
col.min = as.numeric(col.min);
col.max = as.numeric(col.max);
x = (pmax(val.range[1], pmin(val.range[2], x))-val.range[1])/diff(val.range);
col.min = pmax(0, pmin(1, col.min))
col.max = pmax(0, pmin(1, col.max))
if (invert)
{
tmp = col.max
col.max = col.min
col.min = tmp
}
nna = !is.na(x);
out = rep(na.col, length(x))
out[nna] = rgb((col.max[1]-col.min[1])*x[nna] + col.min[1],
(col.max[2]-col.min[2])*x[nna] + col.min[2],
(col.max[3]-col.min[3])*x[nna] + col.min[3])
return(out)
}
#' @name brewer.master
#' @title brewer.master
#' @description
#' Make brewer colors using entire palette
#'
#' Makes a lot of brewer colors using entire brewer palette
#' @param scalar positive integer of number of colors to return
#' @param palette any brewer.pal palette to begin with (default 'Accent')
#'
#' @export
#' @keywords internal
#' @author Marcin Imielinski
brewer.master = function(n, palette = 'Accent')
{
palettes = list(
sequential = c('Blues'=9,'BuGn'=9, 'BuPu'=9, 'GnBu'=9, 'Greens'=9, 'Greys'=9, 'Oranges'=9, 'OrRd'=9, 'PuBu'=9, 'PuBuGn'=9, 'PuRd'=9, 'Purples'=9, 'RdPu'=9, 'Reds'=9, 'YlGn'=9, 'YlGnBu'=9, 'YlOrBr'=9, 'YlOrRd'=9),
diverging = c('BrBG'=11, 'PiYG'=11, 'PRGn'=11, 'PuOr'=11, 'RdBu'=11, 'RdGy'=11, 'RdYlBu'=11, 'RdYlGn'=11, 'Spectral'=11),
qualitative = c('Accent'=8, 'Dark2'=8, 'Paired'=12, 'Pastel1'=8, 'Pastel2'=8, 'Set1'=9, 'Set2'=8, 'Set3'=12)
);
palettes = unlist(palettes);
names(palettes) = gsub('\\w+\\.', '', names(palettes))
if (palette %in% names(palettes))
i = match(palette, names(palettes))
else
i = ((max(c(1, suppressWarnings(as.integer(palette))), na.rm = T)-1) %% length(palettes))+1
col = c();
col.remain = n;
while (col.remain > 0)
{
if (col.remain > palettes[i])
{
next.n = palettes[i]
col.remain = col.remain-next.n;
}
else
{
next.n = col.remain
col.remain = 0;
}
col = c(col, RColorBrewer::brewer.pal(max(next.n, 3), names(palettes[i])))
i = ((i) %% length(palettes))+1
}
col = col[1:n]
return(col)
}
#' @name lighten
#' @title lighten
#' @description
#' lighten
#'
#' lightens / darkens colors by brighness factor f (in -255 .. 255) that will make lighter if > 0 and darker < 0
#' @author Marcin Imielinski
#' @keywords internal
lighten = function(col, f)
{
M = col2rgb(col)
return(apply(matrix(pmax(0, pmin(255, M + f*matrix(rep(1, length(M)), nrow = nrow(M)))), ncol = length(col))/255, 2, function(x) rgb(x[1], x[2], x[3])))
}
#' @name plot.blank
#' @title plot.blank
#' @description
#' Make a blank plot
#'
#' Shortcut for making blank plot with no axes
#' @author Marcin Imielinski
#' @keywords internal
plot.blank = function(xlim = c(0, 1), ylim = c(0,1), xlab = "", ylab = "", axes = F, bg.col = "white", ...)
{
par(bg = bg.col)
plot(0, type = "n", axes = axes, xlab = xlab, ylab = ylab, xlim = xlim, ylim = ylim, ...)
# par(usr = c(xlim, ylim))
}
#' @name draw.triangle
#' @title draw.triangle
#' @description
#' internal draw.triangle function
#'
#' @author Jeremiah Wala
#' @keywords internal
draw.triangle <- function(grl,mdata,y,
ylim.parent=NULL,
windows = NULL,
win.gap = NULL,
#$ m.bg.col = NA, ## DEPRECATED, it is now cmin
sigma = NA, ## if not NA then will blur with a Gaussian filter using a sigma value of this many base pairs
col.min='white',
col.max='red',
gr.colormap = NA,
cmap.min = NA,
cmap.max = NA,
m.sep.lwd = NA,
legend = TRUE,
max.ranges = Inf,
leg.params = list(),
min.gapwidth = 1,
islog = FALSE) {
# if (!is.null(gr.colormap)) {
# if (!is.null(names(gr.colormap)))
# palette.colors = names(gr.colormap)
# }
ylim.subplot = NULL
xlim = c(0, 20000)
if (is.list(y)) {
if (all(c('start', 'end') %in% names(y)))
ylim.subplot = c(y$start[1], y$end[1])
} else {
ylim.subplot = c(y[[1]], y[[2]])
}
if (inherits(grl, "GRangesList")) {
gr <- grl.unlist(grl)
} else {
gr <- grl
}
##assume grl is always non zero
if (is.null(mdata)) {
mdata = matrix(nrow=length(gr), ncol=length(gr), 0)
} else if (nrow(mdata) != length(gr)) {
warning('draw.triangle: square matrix should be same dim as gr')
}
## find covered windows in provided grl
if (is.null(windows)) {
seqlevels(gr) = seqlevels(gr)[seqlevels(gr) %in% as.character(seqnames(gr))]
windows = as(coverage(gr), 'GRanges');
windows = windows[values(windows)$score!=0]
windows = reduce(windows, min.gapwidth = min.gapwidth);
}
if (is.null(win.gap))
win.gap = mean(width(windows))*0.2
## calculate the background
mapped = gr.flatmap(gr, windows, win.gap);
grl.segs = mapped$grl.segs;
window.segs = mapped$window.segs;
winlim = range(c(window.segs$start, window.segs$end + 1)) ## 3/20/16 added +1
#mdata = as.matrix(mdata)[as.numeric(grl.segs$query.id), as.numeric(grl.segs$query.id)]
mdata = mdata[as.numeric(grl.segs$query.id), as.numeric(grl.segs$query.id), drop = FALSE]
#' mimielinski Monday, Jan 28, 2019 12:27:43 AM
if (nrow(grl.segs) > max.ranges)
{
ix = sample(nrow(grl.segs), max.ranges)
grl.segs = grl.segs[ix,]
if (!is.null(mdata))
mdata = mdata[ix, ix, drop = FALSE]
}
if (!is.na(sigma)) ## MARCIN: if blur use spatstat library to blur matrix for n base pairs but making sure we don't bleed across windows
{
uwin = unique(grl.segs$window)
win.map = split(1:length(grl.segs$window), grl.segs$window)
uwin.pairs = cbind(rep(uwin, length(uwin)), rep(uwin, each = length(uwin)))
for (p in 1:nrow(uwin.pairs))
{
print(paste('...blurring', p, 'of', nrow(uwin.pairs), 'windows'))
tmp.i = win.map[[uwin.pairs[p, 1]]]
tmp.j = win.map[[uwin.pairs[p, 2]]]
mdata[tmp.i, tmp.j] = as.matrix(spatstat::blur(spatstat::im(mdata[tmp.i, tmp.j], xcol = grl.segs$start[tmp.j], yrow = grl.segs$end[tmp.i]), sigma = sigma))
}
}
if (nrow(mdata) != nrow(grl.segs))
warning('problem after flatmap. Should have trimmed matrix to len(gr)')
if (nrow(grl.segs) > 0)
if (nrow(grl.segs) != nrow(mdata)) {
warning('problemo')
}
## transform into plot coordinates
grl.segs$pos1 = affine.map(grl.segs$pos1, winlim, ylim = xlim) #d 4 down
grl.segs$pos2 = affine.map(grl.segs$pos2 + 1, winlim, ylim = xlim) ## +1 added
window.segs$start = affine.map(window.segs$start, winlim, ylim = xlim)
window.segs$end = affine.map(window.segs$end + 1, winlim, ylim = xlim) ## 3/20/16 added +1 to make [5,6][7,8] adjacent, etc
#####################
## send to plotting device
#####################
## should have already called draw.grl'!
## make the plot
wid = 20000
hgt.subp = abs(diff(ylim.subplot))
asp.ratio <- (par('pin')[2])/(par('pin')[1])
hgt.plot <- abs(diff(ylim.parent))
dlim <- c(0, wid * asp.ratio * (hgt.subp / hgt.plot))
y0 <- dlim[1]
y1 <- dlim[2]
## draw blank background
rect(xlim[1]-diff(xlim)*0.1, ylim.subplot[1], xlim[2], ylim.subplot[2], border = NA, col = par('bg'))
## set the color scale
if (is.na(cmap.min))
cmap.min = quantile(mdata, 0.01)
if(is.na(cmap.max))
cmap.max = quantile(mdata, 0.99)
#cs <- col.scale(seq(cmap.min, cmap.max), val.range=c(cmap.min, cmap.max), col.min=col.min, col.max=col.max)
if (all(is.na(gr.colormap)))
gr.colormap = c("light green", "yellow", "orange", "red")
else if (is.list(gr.colormap))
gr.colormap = unlist(gr.colormap)
cs <- colorRampPalette(gr.colormap)(length(seq(cmap.min, cmap.max, by=(cmap.max-cmap.min)/100)))
## MARCIN: changed m.bg.col to be equal to cmap.min, ie we just draw it as part of the background
## and save drawing a subset (or majority) of pixels (e.g. those valued 0) if they are just
## "background"
## TODO: what if the "background" pixels are not 0 or cmap.min? would be nice to save on drawing
## then too ...
m.bg.col = cs[1]
## draw the background BOXES using the min color
if (nrow(window.segs) > 1) {
bgx = .all.xpairs(window.segs$start, window.segs$end)
out <- diamond(bgx[,1], bgx[,2], bgx[,3], bgx[,4], y0, y1)
## affine map to local coordinates
out$y[!is.na(out$y)] <- affine.map(out$y[!is.na(out$y)], xlim=dlim, ylim=ylim.subplot)
if (m.sep.lwd > 0)
polygon(out$x, out$y, col=m.bg.col, lwd=m.sep.lwd)
else
polygon(out$x, out$y, col=m.bg.col, border=NA)
}
## draw the background triangles
out = triangle(x1=window.segs$start, x2=window.segs$end, y=y0, y0=y0, y1=y1)
out$y[!is.na(out$y)] <- affine.map(out$y[!is.na(out$y)], xlim=dlim, ylim=ylim.subplot)
#out$y[seq(from=2, to=length(out$y), by=4)] <- affine.map(out$y[seq(from=2, to=length(out$y), by=4)], xlim=dlim, ylim=ylim.subplot)
if (m.sep.lwd > 0)
polygon(out$x, out$y, col=m.bg.col, lwd=m.sep.lwd)
else
polygon(out$x, out$y, col=m.bg.col, border=NA)
if (nrow(grl.segs) == 0)
return(window.segs)
################
## plot the data
################
## set the color scale
#bgx = .all.xpairs(grl.segs$pos1, grl.segs$pos2)
ij = as.data.table(Matrix::which(mdata>cmap.min, arr.ind = TRUE))
setnames(ij, c('i', 'j'))
setkeyv(ij, c('i', 'j'))
bgx = data.frame(grl.segs$pos1[ij$i], grl.segs$pos2[ij$i],
grl.segs$pos1[ij$j], grl.segs$pos2[ij$j],
ij$i, ij$j)
if (nrow(bgx)>0)
{
col = mdata[cbind(bgx[,5], bgx[,6])]
out <- diamond(bgx[,1], bgx[,2], bgx[,3], bgx[,4], y0, y1, col)
## affine map to local coorinates
out$y[!is.na(out$y)] <- affine.map(out$y[!is.na(out$y)], xlim=dlim, ylim=ylim.subplot)
ix.min <- out$col < cmap.min
ix.max <- out$col > cmap.max
out$col[ix.min | ix.max] <- NA
cr <- rep('black', length(out$col))
cr[!is.na(out$col)] <- cs[floor(affine.map(out$col[!is.na(out$col)], xlim = c(cmap.min, cmap.max), ylim=c(1,100)))]
cr[ix.min] <- cs[1] ## deal with cmap.min and cmap.max MARCIN:
cr[ix.max] <- cs[length(cs)] ## deal with cmap.min and cmap.max
# cr[ix.min] <- 'white' ## deal with cmap.min and cmap.max MARCIN:
# cr[ix.max] <- 'black' ## deal with cmap.min and cmap.max
polygon(out$x, out$y, col=cr, border=NA)
}
## plot the triangles part
col = Matrix::diag(mdata)
out.t = triangle(grl.segs$pos1, grl.segs$pos2, y0, y0, y1, col=col)
out.t$y[!is.na(out.t$y)] <-
affine.map(out.t$y[!is.na(out.t$y)], xlim=dlim, ylim=ylim.subplot)
ix.min <- out.t$col < cmap.min
ix.max <- out.t$col > cmap.max
out.t$col[ix.min | ix.max] <- NA
cr <- rep('black', length(out.t$col))
cr[!is.na(out.t$col)] <- cs[floor(affine.map(out.t$col[!is.na(out.t$col)], xlim = c(cmap.min, cmap.max), ylim=c(1,100)))]
cr[ix.min] <- cs[1] ## deal with cmap.min and cmap.max MARCIN:
cr[ix.max] <- cs[length(cs)] ## deal with cmap.min and cmap.max
# cr[ix.min] <- 'white' ## deal with cmap.min and cmap.max
# cr[ix.max] <- 'black' ## deal with cmap.min and cmap.max
#cr <- cs[ceiling(out.t$col) - cmap.min + 1]
polygon(out.t$x, out.t$y, col=cr, border=NA)
legend.cex = leg.params$cex
if (is.null(legend.cex))
legend.cex = 1
if (is.null(legend))
legend = TRUE
if (is.na(legend))
legend = TRUE
if (legend)
{
## plot the legend
if (!islog)
txt = format(c(cmap.min, 0.5*(cmap.max-cmap.min) + cmap.min, cmap.max), digits=1)
else
txt = format(c(10^cmap.min, 10^(0.3333*((cmap.max-cmap.min) + cmap.min)), 10^(0.66667*((cmap.max-cmap.min) + cmap.min)), 10^(cmap.max)), digits=1)
color.bar(lut = cs, xpos=0, ypos=ylim.subplot[1] + diff(ylim.subplot)*0.3, width=500, height=diff(ylim.subplot)*0.3,
text=txt, cex=legend.cex)
}
return(window.segs)
}
#' @name clip_polys
#' @title clip_polys
#' @description
#' Clip polygons
#' @author Jeremiah Wala
#' @keywords internal
clip_polys <- function(dt, y0, y1) {
## fix global def NOTE
y2 <- y3 <- y4 <- y5 <- y6 <- NULL
left <- x3.tmp <- x3 <- y3.tmp <- lean.sign <- NULL
## leave early if not necessary
miny = min(dt[, c(y1, y2, y3, y4, y5, y6)], na.rm = TRUE)
maxy = max(dt[, c(y1, y2, y3, y4, y5, y6)], na.rm = TRUE)
##print(paste('MinY:', miny, 'MaxY:', maxy, "Y0:", y0, 'Y1:', y1))
if (y0 <= miny && y1 >= maxy)
return(dt)
dt[, left := y3 > y6]
dt[, x3.tmp := x3]
dt[, y3.tmp := y3]
dt$x3[dt$left] <- dt$x6[dt$left]
dt$y3[dt$left] <- dt$y6[dt$left]
dt$x6[dt$left] <- dt$x3.tmp[dt$left]
dt$y6[dt$left] <- dt$y3.tmp[dt$left]
dt[, lean.sign := ifelse(left, -1, 1)]
## remove stuff all the way out
yc1 = y1
yc0 = y0
dt <- subset(dt, y1 < yc1 & y4 > yc0)
#############
# deal with the top clip
#############
## clip to only bottom corner
ix <- dt$y3 >= y1
dt$y3[ix] <- dt$y6[ix] <- y1
dt$x3[ix] <- dt$x1[ix] - dt$lean.sign[ix] * (y1 - dt$y1[ix])
dt$x4[ix] <- dt$x5[ix] <- dt$x3[ix]
dt$y4[ix] <- dt$y5[ix] <- dt$y3[ix]
dt$x6[ix] <- dt$x1[ix] + dt$lean.sign[ix] * (y1 - dt$y1[ix])
## clip off 4,5,6
ix <- dt$y6 >= y1 & !ix
dt$y4[ix] <- dt$y5[ix] <- dt$y6[ix] <- y1
dt$x4[ix] <- dt$x5[ix] <- dt$x3[ix] + dt$lean.sign[ix] * (y1 - dt$y3[ix])
dt$x6[ix] <- dt$x1[ix] + dt$lean.sign[ix] * (y1 - dt$y1[ix])
## clip off 4,5
ix <- dt$y4 > y1 & !ix
dt$y4[ix] <- dt$y5[ix] <- y1
dt$x4[ix] <- dt$x3[ix] + dt$lean.sign[ix] * (y1 - dt$y3[ix])
dt$x5[ix] <- dt$x6[ix] - dt$lean.sign[ix] * (y1 - dt$y6[ix])
##################
# deal with the bottom clip
##################
## clip to only top corner
ix <- dt$y6 <= y0
dt$x3[ix] <- dt$x3[ix] + dt$lean.sign[ix] * (y0 - dt$y3[ix])
dt$x6[ix] <- dt$x6[ix] - dt$lean.sign[ix] * (y0 - dt$y6[ix])
dt$y3[ix] <- dt$y6[ix] <- y0
dt$y1[ix] <- dt$y2[ix] <- dt$y3[ix]
dt$x1[ix] <- dt$x2[ix] <- dt$x3[ix]
## clip off 1,2,3
ix <- dt$y3 < y0 & !ix
dt$x3[ix] <- dt$x3[ix] + dt$lean.sign[ix] * (y0 - dt$y3[ix])
dt$x1[ix] <- dt$x1[ix] + dt$lean.sign[ix] * (y0 - dt$y1[ix])
dt$y1[ix] <- dt$y3[ix] <- y0
dt$y2[ix] <- dt$y1[ix]
dt$x2[ix] <- dt$x1[ix]
## clip off 1,2
ix <- dt$y1 < y0 & !ix
dt$x1[ix] <- dt$x1[ix] - dt$lean.sign[ix] * (y0 - dt$y1[ix])
dt$x2[ix] <- dt$x2[ix] + dt$lean.sign[ix] * (y0 - dt$y2[ix])
dt$y1[ix] <- dt$y2[ix] <- y0
return(dt)
}
#' @name diamond
#' @title diamond
#' @description
#' internal draw.triangle function
#' @author Jeremiah Wala
#' @keywords internal
diamond <- function (x11, x12, x21, x22, y0, y1, col=NULL) {
i1 = i2 = .geti(x12, x21)
i3 = .geti(x11, x21)
i4 = i5 = .geti(x11, x22)
i6 = .geti(x12, x22)
## dummy
if (is.null(col))
col <- i1$x
## setup the data table
dt <- data.table::data.table(x1=i1$x, x2=i2$x, x3=i3$x, x4=i4$x, x5=i5$x, x6=i6$x,
y1=i1$y, y2=i2$y, y3=i3$y, y4=i4$y, y5=i5$y, y6=i6$y, col=col)
dt <- clip_polys(dt, y0, y1)
iN <- rep(NA, nrow(dt))
out.x <- as.numeric(t(matrix(c(dt$x1, dt$x2, dt$x3, dt$x4, dt$x5, dt$x6, iN), nrow=nrow(dt))))
out.y <- as.numeric(t(matrix(c(dt$y1, dt$y2, dt$y3, dt$y4, dt$y5, dt$y6, iN), nrow=nrow(dt))))
return(list(x=out.x, y=out.y, col=dt$col))
}
#' @name triangle
#' @title triangle
#' @description
#' internal draw.triangle function
#' @author Marcin Imielinski
#' @keywords internal
triangle <- function(x1, x2, y, y0, y1, col=NULL) {
if (length(y) == 1)
y = rep(y, length(x1))
else if (length(y) != length(x1))
warning('triangle: expecnting length(y) == length(x1) or length(y) == 1')
if (length(x1) != length(x2))
warning('triangle: expecting length(x1) == length(x2)')
i1 = i2 = i3 = list(x=x1, y=y)
i4 = i5 = .geti(x1, x2)
i6 = list(x=x2, y=y)
## dummy
if (is.null(col))
col <- i1$x
dt <- data.table::data.table(x1=i1$x, x2=i2$x, x3=i3$x, x4=i4$x, x5=i5$x, x6=i6$x,
y1=i1$y, y2=i2$y, y3=i3$y, y4=i4$y, y5=i5$y, y6=i6$y, col=col)
dt <- clip_polys(dt, y0, y1)
iN <- rep(NA, nrow(dt))
out.x <- as.numeric(t(matrix(c(dt$x1, dt$x2, dt$x3, dt$x4, dt$x5, dt$x6, iN), nrow=nrow(dt))))
out.y <- as.numeric(t(matrix(c(dt$y1, dt$y2, dt$y3, dt$y4, dt$y5, dt$y6, iN), nrow=nrow(dt))))
return(list(x=out.x, y=out.y, col=col))
}
# @description
# internal draw.triangle function
# @keywords internal
.geti <- function(x1, x2) {
if (length(x1) != length(x2))
stop('x1 and x2 need to be same length')
dt = data.frame(x1, x2)
x = rowMeans(dt)
dt$x1 = (-1) * dt$x1
y = abs(rowMeans(dt)) # / slope
list(x=x, y=y)
}
#' @keywords internal
.getpoints <- function (b1, b2) {
i1 = .geti(b1[1], b2[1])
i2 = .geti(b1[1], b2[2])
i3 = .geti(b1[2], b2[2])
i4 = .geti(b1[2], b2[1])
return(list(x=c(i1[1], i2[1], i3[1], i4[1]), y=c(i1[2], i2[2], i3[2], i4[2])))
}
# internal draw.triangle function
# @keywords internal
.all.xpairs <- function(start, end) {
if (length(start) != length(end)) {
warning('.all.diamonds: lengths of input need to be same. Reducing')
start <- start[min(length(start), length(end))]
end <- end[min(length(start), length(end))]
}
if (length(start) == 1) {
return(matrix())
}
## num.boxes = sum(seq(1, length(start)-1))
## sr = seq(1, length(start)-1)
## rsr = rev(sr)
## ir = rep(sr, rsr)
## jr = unlist(lapply(seq(2,length(start)), function(x) seq(x,length(start))))
## bgx = matrix(c(start[ir], end[ir], start[jr], end[jr], ir, jr), nrow=num.boxes, ncol=6)
## MARCIN speed up, but likely unnecessary
ij = as.data.table(expand.grid(i = 1:length(start), j = 1:length(end)))[i<j, ]
setkeyv(ij, c('i', 'j'))
bgx = matrix(c(start[ij$i], end[ij$i], start[ij$j], end[ij$j], ij$i, ij$j), nrow = nrow(ij), ncol = 6)
return(bgx)
}
#' Function to plot color bar
#' http://stackoverflow.com/questions/9314658/colorbar-from-custom-colorramppalette
#' @keywords internal
color.bar <- function(lut, nticks=11, ticks=seq(min, max, len=nticks), title='',
xpos, ypos, width, height, text=c("min", 'max'), cex=1) {
scale = (length(lut)-1)/(height)
#scale = 1
#dev.new(width=1.75, height=5)
#plot(c(0,10), c(min,max), type='n', bty='n', xaxt='n', xlab='', yaxt='n', ylab='', main=title)
#axis(2, ticks, las=1)
rect(xpos, ypos, xpos+width, ypos + height, lwd=2)
#text(xpos+width, ypos, text[1])
#text(xpos+width, ypos+height, text[2])
for (i in 1:(length(lut)-1)) {
y = (i-1)/scale + ypos
rect(xpos, y, width, y+1/scale, col=lut[i], border=NA)
}
## add the text
for (i in 1:length(text)) {
y = ypos + (i-1)/(length(text)-1) * height
segments(xpos+width, y, xpos+width*1.3, y)
text(xpos+width*2.5, y, text[i], cex=cex)
}
}
#' utility function to load data files into variables
#' @keywords internal
readData = function(..., environment = NULL)
{
my.env = new.env()
data(..., envir = my.env);
return(as.list(my.env))
}
#' gr.flatmap
#'
#' Takes \code{GRanges} pile and maps onto a flattened coordinate system defined by windows (GRanges object)
#' a provided "gap" (in sequence units). If squeeze == TRUE then will additionally squeeze ranges into xlim.
#'
#' output is list with two fields corresponding to data frames:
#' $grl.segs = data frame of input gr's "lifted" onto new flattened coordinate space (NOTE: nrow of this not necessarily equal to length(gr))
#' $window.segs = the coordinates of input windows in the new flattened (and squeezed) space
#' @param gr \code{GRanges} pile to flatten
#' @param windows \code{GRanges} pile of windows defining the coordinate system
#' @param gap TODO
#' @param strand.agnostic TODO
#' @param squeeze TODO
#' @param xlim TODO
#' @return TODO
#' @name gr.flatmap
#' @keywords internal
gr.flatmap = function(gr, windows, gap = 0, strand.agnostic = TRUE, squeeze = FALSE, xlim = c(0, 1))
{
if (strand.agnostic)
GenomicRanges::strand(windows) = "*"
## now flatten "window" coordinates, so we first map gr to windows
## (replicating some gr if necessary)
# h = findOverlaps(gr, windows)
h = gr.findoverlaps(gr, windows);
window.segs = gr.flatten(windows, gap = gap)
grl.segs = BiocGenerics::as.data.frame(gr);
grl.segs = grl.segs[values(h)$query.id, ];
grl.segs$query.id = values(h)$query.id;
grl.segs$window = values(h)$subject.id
grl.segs$start = start(h);
grl.segs$end = end(h);
grl.segs$pos1 = pmax(window.segs[values(h)$subject.id, ]$start,
window.segs[values(h)$subject.id, ]$start + grl.segs$start - start(windows)[values(h)$subject.id])
grl.segs$pos2 = pmin(window.segs[values(h)$subject.id, ]$end,
window.segs[values(h)$subject.id, ]$start + grl.segs$end - start(windows)[values(h)$subject.id])
grl.segs$chr = grl.segs$seqnames
if (squeeze)
{
min.win = min(window.segs$start)
max.win = max(window.segs$end)
grl.segs$pos1 = affine.map(grl.segs$pos1, xlim = c(min.win, max.win), ylim = xlim)
grl.segs$pos2 = affine.map(grl.segs$pos2, xlim = c(min.win, max.win), ylim = xlim)
window.segs$start = affine.map(window.segs$start, xlim = c(min.win, max.win), ylim = xlim)
window.segs$end = affine.map(window.segs$end, xlim = c(min.win, max.win), ylim = xlim)
}
return(list(grl.segs = grl.segs, window.segs = window.segs))
}
gr.stripstrand = function(gr)
{
GenomicRanges::strand(gr) = "*"
return(gr)
}
format_windows <- function(windows, .Object) {
if (is(windows, 'character'))
windows = BiocGenerics::unlist(parse.grl(windows, seqlengths(seqinfo(.Object))))
if (length(windows)==0)
stop('Empty windows provided or error parsing windows - may want to check intervals against DEFAULT_GENOME environment variable or output of hg_seqlengths()')
if (is(windows, 'Seqinfo'))
windows = si2gr(windows)
if (is(windows, 'GRangesList'))
windows = BiocGenerics::unlist(windows)
windows = windows[width(windows)>0] ## otherwise will get non-matching below
if (is.null(windows$col))
windows$col = 'gray98'
# if (is.list(windows))
# windows = do.call('GRangesList', windows)
## collapse and match metadata back to original
tmp = reduce(gr.stripstrand(windows))
ix = gr.match(tmp, windows)
values(tmp) = values(windows)[ix, , drop = FALSE]
tmp = tmp[order(ix)] ## try to respect initial order
windows = tmp
## if (!inherits(windows, 'GRangesList')) ## GRangesList windows deprecated
## windows = GenomicRanges::GRangesList(windows)
if (sum(as.numeric(width(grl.unlist(windows))))==0)
{
if (length(seqinfo(.Object))) {
windows = si2gr(seqinfo(.Object))
} else {
warning("no windows provided and no seqinfo. Drawing blank plot")
return(GRanges())
}
}
return (grl.unlist(windows))
}
prep_defaults_for_plotting <- function(.Object) {
# if (is.null(.Object@formatting$triangle))
# .Object@formatting$triangle = NA
#
# if (any(ix <- is.na(.Object@formatting$triangle)))
# .Object@formatting$triangle[ix] = FALSE
#
if (is.null(formatting(.Object)$source.file.chrsub))
formatting(.Object)$source.file.chrsub = TRUE
# layout legends if colorfield or colormap is NA and legends have no xpos set
# leg.ix = which(.Object@formatting$legend & (!is.na(.Object@formatting$gr.colorfield) | !sapply(.Object@colormap, is.null)))
# if (is.null(.Object@formatting$legend.xpos))
# .Object@formatting$legend.xpos = NA
#
# if (is.null(.Object@formatting$legend.xjust))
# .Object@formatting$legend.xjust = NA
#
# if (any(is.na(.Object@formatting$legend.xpos[leg.ix])))
# {
# .Object@formatting$legend.xpos[leg.ix] = seq(0.1, 1, length.out = length(leg.ix))
# .Object@formatting$legend.xjust[leg.ix] = ifelse(.Object@formatting$legend.xpos[leg.ix]<0.5, 0, 1)
# }
# layout y coordinates
sumh = sum(formatting(.Object)$height + formatting(.Object)$ygap)
formatting(.Object)$height = formatting(.Object)$height/sumh
formatting(.Object)$ygap = formatting(.Object)$ygap/sumh
# formatting(.Object)$ywid = formatting(.Object)$ywid/sumh
if (is.null(formatting(.Object)$max.ranges))
formatting(.Object)$max.ranges = NA
return(.Object)
}
extract_data_from_tmp_dat <- function(.Object, j, this.windows) {
tmp.dat = dat(.Object)[[j]]
# if (inherits(tmp.dat, 'RleList'))
# {
# tmp.score = rle.query(tmp.dat, this.windows)
# tmp.dat = gUtils::gr.dice(this.windows)
# tmp.dat$score = as.numeric(tmp.score)
#
# if (is.na(formatting(.Object)$y0[j]))
# formatting(.Object)$y0[j] = min(tmp.dat$score, na.rm = T)
#
# if (is.na(formatting(.Object)$y1[j]))
# formatting(.Object)$y1[j] = max(tmp.dat$score, na.rm = T)
#
# formatting(.Object)$y.field[j] = 'score'
# pre.filtered = TRUE
# }
## else if
if (is.character(tmp.dat))
{
if (file.exists(tmp.dat))
{
bed.style = FALSE
if (grepl('(\\.bw)|(\\.bigwig)', tmp.dat, ignore.case = TRUE))
{
f = rtracklayer::BigWigFile(normalizePath(tmp.dat))
si = tryCatch(seqinfo(f), error = function(tmp.dat) NULL)
}
else if (grepl('\\.wig', tmp.dat, ignore.case = TRUE))
{
f = rtracklayer::WIGFile(normalizePath(tmp.dat))
si = tryCatch(seqinfo(f), error = function(tmp.dat) NULL)
}
else if (grepl('\\.bed', tmp.dat, ignore.case = TRUE))
{
f = rtracklayer::BEDFile(normalizePath(tmp.dat))
# si = tryCatch(seqinfo(f), error = function(tmp.dat) NULL)
bed.style = TRUE
}
else if (grepl('\\.gff', tmp.dat, ignore.case = TRUE))
{
f = rtracklayer::GFFFile(normalizePath(tmp.dat))
si = tryCatch(seqinfo(f), error = function(tmp.dat) NULL)
}
else if (grepl('\\.2bit', tmp.dat, ignore.case = TRUE))
{
f = rtracklayer::TwoBitFile(normalizePath(tmp.dat))
si = tryCatch(seqinfo(f), error = function(tmp.dat) NULL)
}
else if (grepl('\\.bedgraph', tmp.dat, ignore.case = TRUE))
{
f = rtracklayer::BEDGraphFile(normalizePath(tmp.dat))
# si = tryCatch(seqinfo(f), error = function(tmp.dat) NULL)
bed.style = TRUE
}
if (!bed.style) ## bed.style file objects do not have seqinfo option
{
si = tryCatch(seqinfo(f), error = function(x) '')
pre.filtered = T
if (!is.character(si))
{
if (formatting(.Object)[j, 'source.file.chrsub'])
{
sel = gr.fix(gr.chr(this.windows), si, drop = TRUE)
}
else
sel = gr.fix(this.windows, si, drop = TRUE)
if (length(this.windows)>0 & length(sel)==0)
warning('zero ranges intersecting UCSC file sequences selection perhaps source.file.chrsub parameter needs to be fixed?')
tmp.dat = rtracklayer::import(f, selection = sel, asRangedData = FALSE)
if (is.na(formatting(.Object)$y.field[j]))
formatting(.Object)$y.field[j] = 'score'
}
}
else
{
tmp.dat = rtracklayer::import(f, asRangedData = FALSE)
if (is.na(formatting(.Object)$y.field[j]))
formatting(.Object)$y.field[j] = 'score'
}
}
}
else if (is(tmp.dat, 'ffTrack'))
{
tmp.dat = tmp.dat[this.windows, gr = TRUE]
formatting(.Object)$y.field[j] = 'score'
pre.filtered = T
}
else
{
if (formatting(.Object)[j, 'source.file.chrsub'])
tmp.dat = gUtils::gr.sub(tmp.dat, 'chr', '')
}
if (is.character(tmp.dat)) ## file was not found
{
warnings('Track bigwig file not found')
tmp.dat = GRanges()
}
if (formatting(.Object)[j, 'source.file.chrsub'])
{
tmp.dat = gUtils::gr.sub(tmp.dat, 'chr', '')
this.windows= gUtils::gr.sub(this.windows, 'chr', '')
}
return(list(o=.Object, t=tmp.dat, w = this.windows))
}
enforce_max_ranges <- function(.Object, pre.filtered, j, tmp.dat, this.windows) {
## enforce max.ranges (everything should be a GRanges at this point)
if (!pre.filtered & nrow(edgs(.Object)[[j]])==0 & length(.Object@vars[[j]])==0) ## assume any track associated with edges is pre-filtered
{
if (length(tmp.dat)>formatting(.Object)$max.ranges[j])
if (inherits(tmp.dat, 'GRangesList'))
{
vals = values(tmp.dat)
nm = names(tmp.dat)
tmp2 = grl.unlist(tmp.dat)
tmp2 = tmp2[gUtils::gr.in(tmp2, this.windows)]
##tmp2 = tmp2[tmp2 %over% this.windows]
tmp.dat = GenomicRanges::split(tmp2, tmp2$grl.ix)
values(tmp.dat) = vals[as.numeric(names(tmp.dat)), ]
names(tmp.dat) = nm[as.numeric(names(tmp.dat))]
}
else if (length(tmp.dat)>formatting(.Object)$max.ranges[j])
{
tmp.dat = tmp.dat[gUtils::gr.in(tmp.dat, this.windows)]
##tmp.dat = tmp.dat[tmp.dat %over% this.windows]
pre.filtered = TRUE
}
}
if (length(tmp.dat)>formatting(.Object)$max.ranges[j] &
nrow(edgs(.Object)[[j]])==0 & !.Object@formatting$triangle[j]) ## don't sample if there are edges or triangle
{
tmp.dat = sample(tmp.dat, ceiling(formatting(.Object)$max.ranges[j]))
}
return(list(p=pre.filtered, t=tmp.dat))
}
smooth_yfield <- function(.Object, j, tmp.dat) {
tmp = GenomicRanges::coverage(tmp.dat, weight = values(tmp.dat)[, formatting(.Object)$y.field[j]])
## calculate footprint of input ranges so we only smooth that region + $smooth parameter window around it
## so as to save time not smoothing entire genome worth of 0s (not clear why feature not already included in S4 vectors runmean)
tmpfp = as(GenomicRanges::coverage(tmp.dat), 'GRanges')
tmpfp = gr2dt(reduce(tmpfp[tmpfp$score!=0])+formatting(.Object)$smooth[j]) ##
tmpfp[, seqnames := as.character(seqnames)]
tmpfp[, end := pmin(lengths(tmp)[as.character(seqnames)], end)]
tmpfp[, start := pmax(start, 1)]
for (i in 1:nrow(tmpfp))
{
tmp[[tmpfp[i, seqnames]]][tmpfp[i, as.vector(start:end)]] =
S4Vectors::runmean(tmp[[tmpfp[i, seqnames]]][tmpfp[i, as.vector(start:end)]],
k = floor(formatting(.Object)$smooth[j]/2)*2+1, endrule = 'drop', na.rm = TRUE)
}
if (!is.na(formatting(.Object)$round[j]))
tmp = round(tmp, formatting(.Object)$round[j])
tmp = as(tmp, 'GRanges')
tmp = tmp[gUtils::gr.in(tmp, tmp.dat)]
##tmp = tmp[tmp %over% tmp.dat]
tmp.val = tmp$score
values(tmp) = values(tmp.dat)[gr.match(tmp, tmp.dat), , drop = F]
values(tmp)[, formatting(.Object)$y.field[j]] = tmp.val
tmp.dat = tmp
}
format_yfield_limits <- function(.Object, j, tmp.dat, pre.filtered, this.windows) {
strand(this.windows) <- rep("*", length(this.windows))
if (!(formatting(.Object)$y.field[j] %in% names(values(tmp.dat))))
stop('y.field missing from input granges')
y0.global = min(values(tmp.dat)[, formatting(.Object)$y.field[j]], na.rm = TRUE)
y1.global = max(values(tmp.dat)[, formatting(.Object)$y.field[j]], na.rm = TRUE)
if (!pre.filtered)
{
if (inherits(tmp.dat, 'GRanges'))
{
tmp.dat.r <- tmp.dat[gUtils::gr.in(tmp.dat, this.windows)]
}
else if (inherits(tmp.dat, 'GRangesList'))
{
nm = names(tmp.dat)
val = values(tmp.dat)
tmp = grl.unlist(tmp.dat)
nmcol = setdiff(names(values(tmp)), names(val))
tmp2 = tmp[gUtils::gr.in(tmp, this.windows)]
tmp.dat.r = split(tmp2[, nmcol], tmp2$grl.ix)
values(tmp.dat.r) = val[as.numeric(names(tmp.dat.r)), , drop = FALSE]
names(tmp.dat.r) = nm[as.numeric(names(tmp.dat.r))]
}
else
{
stop('Unrecognized format inside gTrack')
}
}
else
{
tmp.dat.r = tmp.dat
}
val = values(tmp.dat.r)[, formatting(.Object)$y.field[j]]
# r = range(val[!is.infinite(val)], na.rm = TRUE)
p.quantile = 0.01
if (!is.null(formatting(.Object)$y.quantile[j]) && !is.na(formatting(.Object)$y.quantile[j]))
p.quantile = pmin(pmax(pmin(formatting(.Object)$y.quantile[j], 1-formatting(.Object)$y.quantile[j]), 0), 1)
r = quantile(val[!is.infinite(val)], probs = c(p.quantile, 1-p.quantile), na.rm = TRUE)
if (is.na(diff(r)))
r = c(0, 0)
if (is.na(formatting(.Object)$y0[j])) ## adjust y limits here if not specified
{
formatting(.Object)$y0[j] = pmax(y0.global, r[1] - diff(r)*0.05)
if (diff(r) == 0 & !is.na(formatting(.Object)$y1)[j])
formatting(.Object)$y0[j] =
r[1] - diff(c(r[1], formatting(.Object)$y1[j]))*0.05
}
if (is.na(formatting(.Object)$y1[j])) ## adjust y limits here if not specified
{
formatting(.Object)$y1[j] = pmin(y1.global, r[2] + diff(r)*0.05)
if (diff(r) == 0 & !is.na(formatting(.Object)$y0[j]))
formatting(.Object)$y1[j] =
r[2] + diff(c(formatting(.Object)$y0[j], r[1]))*0.05
}
if (!is.null(formatting(.Object)$y0.bar[j]) && !is.na(formatting(.Object)$y0.bar[j]))
formatting(.Object)$y0[j] = min(y0.global, formatting(.Object)$y0.bar[j])
return(.Object)
}
draw_x_ticks <- function(xaxis.interval = 'auto', windows, mapped, winlim, xlim, ylim,
xaxis.pos = 1,
xaxis.suffix = 'MB',
xaxis.unit = 1,
xaxis.cex.tick = 1,
xaxis.ticklen = 1,
xaxis.round = 3
) {
wid <- sum(as.numeric(width(windows)))
xaxis.nticks = 20 ## default number of ticks
if (is.na(xaxis.unit))
xaxis.unit = 'MB'
if (is.na(xaxis.cex.tick))
xaxis.cex.tick =1
if (is.na(xaxis.ticklen))
xaxis.ticklen = 1
if (is.na(xaxis.round))
xaxis.round = 1
if (is.na(xaxis.pos))
xaxis.pos = 1
if (is.na(xaxis.interval))
xaxis.interval = 'auto'
if (xaxis.interval == 'auto') {
if (wid > 100)
xaxis.interval = 10^(ceiling(log10(wid/xaxis.nticks)))
else
xaxis.interval <- max(floor(wid/xaxis.nticks), 1)
}
# don't let xaxis.interval to be too small .. so tick drawing doesn't get too out of control
xaxis.nticks = 1000 ## max number of ticks
if (!is.na(xaxis.interval))
xaxis.interval = max(xaxis.interval, 10^(ceiling(log10(sum(as.numeric(width(windows)))/xaxis.nticks))))
seq.at.og = lapply(1:length(windows), function(x)
{
out = c(start(windows)[x], seq(ceiling(start(windows)[x]/xaxis.interval),
floor(end(windows)[x]/xaxis.interval))*xaxis.interval, end(windows)[x])
out[out>=start(windows)[x] & out<=end(windows)[x]]
});
winlim[2] <- winlim[2] + 1
seq.at = lapply(1:length(seq.at.og), function(x) affine.map(seq.at.og[[x]]-start(windows)[x]+mapped$window.segs$start[x] + 0.5, ## added 0.5
xlim = winlim, ylim = xlim))
seq.at.og = unlist(seq.at.og)
seq.at = unlist(seq.at);
dup.ix = duplicated(seq.at);
seq.at.og = seq.at.og[!dup.ix]
#seq.at.og <- seq.at.og[-length(seq.at.og)]## jeremiah
seq.at = seq.at[!dup.ix]
# then ticks
seq.at <- seq.at[!is.na(seq.at)]
tick.len = 0.01*diff(ylim)*xaxis.ticklen
y0.tick = rep(xaxis.pos, length(seq.at))
y1.tick = rep(xaxis.pos, length(seq.at)) - tick.len
segments(seq.at[!is.na(seq.at)], y0.tick, seq.at, y1.tick)
# then (tick) text
if (xaxis.unit == 1)
tick.text = prettyNum(paste(seq.at.og, xaxis.suffix), big.mark = ',')
else
tick.text = paste(round(seq.at.og/xaxis.unit, xaxis.round), xaxis.suffix)
##if (xaxis.nticks > 0)
if (!is.na(xaxis.interval) && xaxis.interval > 0)
text(seq.at, y1.tick-tick.len, tick.text,
cex = xaxis.cex.tick*0.6, srt = 90, adj = c(1, 0.5))
}
#' \code{stats::aggregate}, but returns vector
#'
#' @description
#' Same as \code{stats::aggregate} except returns named vector
#' with names as first column of output and values as second
#'
#' Note: there is no need to ever use aggregate or vaggregate, just switch to data.table
#'
#' @param ... arguments to aggregate
#' @return named vector indexed by levels of "by"
#' @author Marcin Imielinski
#' @keywords internal
vaggregate = function(...)
{
out = aggregate(...);
return(structure(out[,ncol(out)], names = do.call(paste, lapply(names(out)[1:(ncol(out)-1)], function(x) out[,x]))))
}
# clip.seg
#
# (data frame seg function)
#
# Clips a pile of segments / ranges to a given region (a single seg) .. ie removing all segments outside of that region
# and then trimming the segments overlapping "region" to the region boundaries
#
clip.seg = function(x, region)
{
x = standardize_segs(x);
region = standardize_segs(region);
if (nrow(region)>1)
stop('Can only region to a single range')
## if either x or region do not specify chrom then we ignore chrom
if (is.null(x$chr))
region$chr = NULL;
if (is.null(region$chr))
x$chr = NULL;
ix = which(seg.on.seg(x, region));
x = x[ix, ]
if (nrow(x)>0)
{
x$pos1 = pmax(region$pos1, x$pos1)
x$pos2 = pmin(region$pos2, x$pos2)
}
return(x)
}
#' standardize_segs
#'
#' (data frame seg function)
#'
#' Takes and returns segs data frame standardized to a single format (ie $chr, $pos1, $pos2)
#'
#' if chr = TRUE will ensure "chr" prefix is added to chromossome(if does not exist)
#' @keywords internal
standardize_segs = function(seg, chr = FALSE)
{
#if (inherits(seg, 'IRangesList'))
# seg = irl2gr(seg);
if (is(seg, 'matrix'))
seg = as.data.frame(seg, stringsAsFactors = FALSE)
if (inherits(seg, 'RangedData') | inherits(seg, 'GRanges') | inherits(seg, 'IRanges'))
{
val = as.data.frame(values(seg));
values(seg) = NULL;
seg = as.data.frame(seg, row.names = NULL); ## returns compressed iranges list
seg$seqnames = as.character(seg$seqnames)
}
else
val = NULL;
field.aliases = list(
ID = c('id', 'patient', 'Sample'),
chr = c('seqnames', 'chrom', 'Chromosome', "contig", "seqnames", "seqname", "space", 'chr', 'Seqnames'),
pos1 = c('start', 'loc.start', 'begin', 'Start', 'start', 'Start.bp', 'Start_position', 'pos', 'pos1', 'left', 's1'),
pos2 = c('end', 'loc.end', 'End', 'end', "stop", 'End.bp', 'End_position', 'pos2', 'right', 'e1'),
strand = c('strand', 'str', 'strand', 'Strand', 'Str')
);
if (is.null(val))
val = seg[, setdiff(names(seg), unlist(field.aliases))]
seg = seg[, intersect(names(seg), unlist(field.aliases))]
for (field in setdiff(names(field.aliases), names(seg)))
if (!(field %in% names(seg)))
names(seg)[names(seg) %in% field.aliases[[field]]] = field;
if (chr)
if (!is.null(seg$chr))
if (!grepl('chr', seg$chr[1]))
seg$chr = paste('chr', seg$chr, sep = "");
if (is.null(seg$pos2))
seg$pos2 = seg$pos1;
missing.fields = setdiff(names(field.aliases), c(names(seg), c('chr', 'ID', 'strand')));
if (length(missing.fields)>0)
warning(sprintf('seg file format problem, missing an alias for the following fields:\n\t%s',
paste(sapply(missing.fields, function(x) paste(x, '(can also be', paste(field.aliases[[x]], collapse = ', '), ')')), collapse = "\n\t")));
if (!is.null(val))
seg = cbind(seg, val)
return(seg)
}
# seg.on.seg
#
# (data frame version of GenomicRanges %in%)
#
# Given two collections of segs segs1 and segs2 (each a data fram of $chr, $start, $end) returns a logical vector of length nrow(segs1)
# that is TRUE in position i if segment i in segs1 overlaps one or more segments in segs2
#
# also works on data frames with the following segment nomenclatures, however data frame should only contain one of these:
# $chrom, $loc.start, $loc.end, $pos1, $pos2
#
# if $chr or $chrom not included then will only consider start and end
#
# e.g. given segments in data frame segs and snps in data frame map snp.in.seg(snps, segs) will flag snps that are in segment
seg.on.seg = function(segs1, segs2, pad = 0, within = F # within = T -> output only T if segs1 range is within at least one segs2 range
)
{
# convert col names to a standard format
segs1 = standardize_segs(segs1);
segs2 = standardize_segs(segs2);
out = vector(mode = "logical", length = nrow(segs1))
# if neither of the sets of segs have chr specified then we act as if all segs are on a single chromosome
if (is.null(segs1$chr) & is.null(segs2$chr))
segs1$chr = segs2$chr = 1;
for (chr in unique(intersect(segs1$chr, segs2$chr)))
{
chrind1 = which(segs1$chr == chr)
chrind2 = which(segs2$chr == chr)
for (i in chrind2)
{
region = segs2[i,, drop = FALSE]
if (within)
out[chrind1] = out[chrind1] | (segs1$pos1[chrind1]>=(region$pos1-pad) & segs1$pos2[chrind1]<=(region$pos2+pad)) # seg1 is inside region
else
out[chrind1] = out[chrind1] | !(segs1$pos2[chrind1]<(region$pos1-pad) | segs1$pos1[chrind1]>(region$pos2+pad))
}
}
return(out)
}
#' Deduplicate character vectors
#
#' Relabels duplicates in a character vector with .1, .2, .3
#' (where "." can be replaced by any user specified suffix)
#' @param x Character vector to deduplicate.
#' @param suffix User defined suffix
#' @return dedupclicated character vector
#' @keywords internal
dedup = function(x, suffix = '.')
{
dup = duplicated(x);
udup = unique(x[dup])
udup.ix = lapply(udup, function(y) which(x==y));
udup.suffices = lapply(udup.ix, function(y) c('', paste(suffix, 2:length(y), sep = '')))
out = x;
out[unlist(udup.ix)] = paste(out[unlist(udup.ix)], unlist(udup.suffices), sep = '');
return(out)
}
get_seqinfo <- function(.Object, seqinfo) {
if (is.na(seqinfo) & is.list(.Object@data)) {
slen = c()
if (is.null(formatting(.Object)$yaxis))
.Object$yaxis = NA
for (i in 1:length(.Object@data))
{
x = .Object@data[[i]]
if (is(x, 'GRanges') || is(x, 'GRangesList'))
{
slen = GenomeInfoDb::seqlengths(x)
if (any(is.na(slen)))
{
if (is(x, 'GRanges'))
slen = GenomeInfoDb::seqlengths(gr.fix(x))
else if (is(x, 'GRangesList'))
slen = GenomeInfoDb::seqlengths(gr.fix(unlist(x)))
}
}
else if (is(x, 'ffTrack'))
{
if (is.null(slen))
slen = GenomeInfoDb::seqlengths(x)
else
slen[seqlevels(x)] = pmax(slen[seqlevels(x)], GenomeInfoDb::seqlengths(x), na.rm = TRUE)
formatting(.Object)[i, 'yaxis'] = TRUE
}
else if (inherits(x, 'RleList'))
{
.Object@data[[i]] = as(x, 'RleList')
formatting(.Object)[i, 'yaxis'] = TRUE
slen = structure(sapply(x, length), names = names(x))
}
else if (is.character(x))
{
if (!file.exists(x))
stop('External file does not exist. Note the file must be a supported UCSC file format or R .rds file with the associated file extension .bw, .wig, .bed., .gff, .2bit, .bedgraph, .rds.')
if (grepl('(\\.bw)|(\\.bigwig)', x, ignore.case = TRUE))
{
f = rtracklayer::BigWigFile(normalizePath(x))
slen = tryCatch(GenomeInfoDb::seqlengths(f), error = function(x) NULL)
formatting(.Object)[i, 'yaxis'] = TRUE
if (is.null(slen))
stop('External file must be a valid and existing .bigwig file')
}
else if (grepl('\\.wig', x, ignore.case = TRUE))
{
f = rtracklayer::WIGFile(normalizePath(x))
slen = tryCatch(GenomeInfoDb::seqlengths(f), error = function(x) NULL)
formatting(.Object)[i, 'yaxis'] = TRUE
if (is.null(slen))
stop('External file must be a valid and existing .wig file')
}
else if (grepl('\\.bed$', x, ignore.case = TRUE)) ## only option is to load
{
f = fread(normalizePath(x))
colnames(f) = c("chrom", "start", "end")
tmp.out = dt2gr(f)
if (is.null(tmp.out))
stop('External file must be a valid and existing .bed file')
else
{
if (.Object@formatting$source.file.chrsub)
.Object@data[[i]] = gUtils::gr.sub(tmp.out, 'chr', '')
else
.Object@data[[i]] = tmp.out
}
slen = GenomeInfoDb::seqlengths(tmp.out)
# slen = tryCatch(seqlengths(f), error = function(x) NULL)
# if (is.null(slen))
# stop('External file must be a valid and existing .bed file')
}
else if (grepl('\\.gff', x, ignore.case = TRUE))
{
f = rtracklayer::GFFFile(normalizePath(x))
slen = tryCatch(GenomeInfoDb::seqlengths(f), error = function(x) NULL)
if (is.null(slen))
stop('External file must be a valid and existing .gff file with a seqlength')
}
else if (grepl('\\.2bit', x, ignore.case = TRUE)) {
f = rtracklayer::TwoBitFile(normalizePath(x))
slen = tryCatch(GenomeInfoDb::seqlengths(f), error = function(x) NULL)
if (is.null(slen))
stop('External file must be a valid and existing .2bit file')
}
else if (grepl('\\.bedgraph', x, ignore.case = TRUE)) ## only option is to load
{
f = rtracklayer::BEDGraphFile(normalizePath(x))
tmp.out = tryCatch(rtracklayer::import(f), error = function(x) NULL)
if (is.null(tmp.out))
stop('External file must be a valid and existing .bedgraph file')
else
{
if (.Object@formatting$source.file.chrsub)
.Object@data[[i]] = gUtils::gr.sub(tmp.out, 'chr', '')
else
.Object@data[[i]] = tmp.out
}
slen = GenomeInfoDb::seqlengths(tmp.out)
}
else if (grepl('\\.rds', x, ignore.case = TRUE)) ## assume this is GRanges or GRangesList
{
tmp.out = tryCatch(readRDS(x), error = function(x) NULL)
if (is.null(tmp.out))
stop('External file must be a valid and existing .rds file containing GRanges or GRangesList')
else if (!inherits(tmp.out, 'GRanges') & !inherits(tmp.out, 'GRangesList'))
stop('External file must be a valid and existing .rds file containing GRanges or GRangesList')
else
{
if (.Object@formatting$source.file.chrsub)
.Object@data[[i]] = gr.sub(tmp.out, 'chr', '')
else
.Object@data[[i]] = tmp.out
}
slen = GenomeInfoDb::seqlengths(tmp.out)
}
}
else
if (is.null(slen))
slen= structure(sapply(x, length), names = names(x))
else
slen[names(x)] = pmax(slen[x], sapply(x, length), na.rm = TRUE)
}
if (any(.Object@formatting$chr.sub))
names(slen) = gsub('chr', '', names(slen))
if (any(duplicated(names(slen))))
slen = data.table(len = slen, nm = names(slen))[, list(len = max(len)), by = nm][, structure(len, names = nm)]
.Object@seqinfo = Seqinfo(seqnames = names(slen), seqlengths = slen);
} else {
if (is.null(seqinfo))
.Object@seqinfo = Seqinfo()
else if (is.list(seqinfo) & all(sapply(seqinfo, inherits, 'Seqinfo')) & length(seqinfo) == length(.Object))
.Object@seqinfo = seqinfo
else
{
if (!inherits(seqinfo, 'Seqinfo'))
seqinfo = GenomeInfoDb::seqinfo(seqinfo);
.Object@seqinfo = seqinfo;
.Object@data = lapply(.Object@data, gr.fix, seqinfo)
}
}
return(.Object)
}
## #' @importFrom data.table as.data.table
## #' @importFrom gUtils dt2gr
## #' @name track.straw
## #' @title track.straw
## #' @description
## #' queries .hic object via straw API https://github.com/theaidenlab/straw/tree/master/R
## #' to extract all of the data in length n query gr (note will do n choose 2 queries since
## #' straw only supports pairwise queries)
## #'
## #'
## #' @keywords straw
## #' @param hic path to .hic file
## #' @param gr granges to query
## #' @param norm string specifying normalization to apply ("KR" (default), "NONE", VC")
## #' @param res resolution to query (default 10kb)
## #' @param type "BP" vs "FRAG"
## #' @param mc.cores parallelization to apply for query (useful when length(gr)>2)
## #' @return gTrack gt of output with mdata(gt)[[1]] is the matrix of contact info and dat(gt)[[1]] is the granges
## #' @export
## #' @author Marcin Imielinski
## track.straw = function(hic, gr, norm = "KR", type = 'BP', res = 1e4, mc.cores = 1, colormap = c('white', 'red', 'black'), ...)
## {
## gr = reduce(gr.stripstrand(gr[, c()]))
## grs = as.data.table(gr)[, paste(seqnames, start, end, sep = ":")]
## n = length(gr)
## grs.pairs = chr2.pairs = chr1.pairs = NULL
## grs.singles = paste(grs, grs)
## chr1.singles = chr2.singles = as.character(seqnames(gr))
## if (n>1)
## {
## pairs = t(combn(n, 2)) ##pairs
## grs.pairs = paste(grs[pairs[,1]], grs[pairs[,2]]) ## combinations
## #' Tuesday, Nov 07, 2017 04:43:52 PM - Julie: adding as.character()
## chr1.pairs = as.character(seqnames(gr)[pairs[,1]])
## chr2.pairs = as.character(seqnames(gr)[pairs[,2]])
## }
## str = paste(norm, hic, c(grs.singles, grs.pairs), type, as.integer(res))
## chr1 = as.character(c(chr1.singles, chr1.pairs))
## chr2 = as.character(c(chr2.singles, chr2.pairs))
## out = rbindlist(mcmapply(str = str, chr1 = chr1, chr2 = chr2,
## FUN = function(str, chr1, chr2)
## {
## dt = as.data.table(.Call("_gTrack_straw_R", str))
## dt[, ":="(chr1 = chr1, chr2 = chr2)]
## return(dt)
## }, SIMPLIFY = FALSE, mc.cores = mc.cores))
## out = out[!is.na(counts), ]
## out[, x2 := x+res-1]
## out[, y2 := y+res-1]
## out[, str1 := paste0(chr1, ':', x, '-', x2)]
## out[, str2:= paste0(chr2, ':', y, '-', y2)]
## gr.out = unique(dt2gr(rbind(out[, .(seqnames = chr1, start = x, end = x2)],
## out[, .(seqnames = chr2, start = y, end = y2)])))
## gr.out$grs = gr.string(gr.out)
## out[, i1 := match(str1, gr.out$grs)]
## out[, j1 := match(str2, gr.out$grs)]
## out[, i := pmin(i1, j1)]
## out[, j := pmax(i1, j1)]
## mdata = Matrix::sparseMatrix(out$i, out$j, x = out$counts, dims = c(length(gr.out), length(gr.out)),symmetric = TRUE)
## gt = gTrack(gr.out, mdata = mdata, colormap = colormap, ...)
## return(gt)
## }
## convert input data into a list of length 'len' of type FUN
listify <- function(x, FUN, len = 1) {
if (is.null(x))
return(rep(list(FUN()), len))
if (!is.list(x))
return(list(x))
return(x)
}
alpha = function(col, alpha)
{
col.rgb = col2rgb(col)
out = rgb(red = col.rgb['red', ]/255, green = col.rgb['green', ]/255, blue = col.rgb['blue', ]/255, alpha = alpha)
names(out) = names(col)
return(out)
}
##########
##########
##########
setGeneric('with')
setMethod("with", signature(data = "gTrack"), NULL)
setMethod("with", signature(data = "gTrack"), function(data, expr) {
df = as.data.frame(formatting(data))
eval(substitute(expr, parent.frame()), df, parent.frame())
p})
setGeneric('within')
setMethod("within", signature(data = "gTrack"), NULL)
setMethod("within", signature(data = "gTrack"), function(data, expr) {
e = list2env(as.list(formatting(data)))
eval(substitute(expr, parent.frame()), e, parent.frame())
formatting(data) = as.data.frame(as.list(e))[, c(colnames(formatting(data))),drop = FALSE]
return(data)
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.