R/tree-force-layout.R
In farrellja/URD: URD
Defines functions
treeForceStretchCoords
treeForceStretchCoords
treeForceTranslateCoords
treeForceRotateCoords
translateCoords3d
rotateCoords3d
plotTreeForceStore3DView
fdlDensity
treeForcePositionLabels
plotTreeForceDual
plotTreeForce
plotTreeForce2D
treeForceDirectedLayout
Documented in
fdlDensity
plotTreeForce
plotTreeForce2D
plotTreeForceDual
plotTreeForceStore3DView
rotateCoords3d
translateCoords3d
treeForceDirectedLayout
treeForcePositionLabels
treeForceRotateCoords
treeForceStretchCoords
treeForceTranslateCoords
#' Generate force-directed layout using tip-walk data.
#'
#' This function generates a k-nearest neighbor network for cells, based on their visitation
#' frequency by the biased random walks from different tips (and optionally pseudotime), and
#' uses it as input into a force-directed layout (powered by igraph). The force-directed
#' layout is generated in 2 dimensions, and pseudotime is used as a third dimension.
#'
#' Several settings adjust the k-nearest neighbor network to assist in the layout. Outlier
#' cells (that are only closely connected to a single other cells) or outlier edges (with
#' unusually long distances) can be eliminated, though these parameters are disabled by
#' default. However, the dendrogram structure recovered by URD is used by default to
#' refine the k-nearest neighbor network, breaking links between cells that are distant
#' in the dendrogram. This emphasizes producing a parseable layout at the expense of
#' representing rare transitions in the data.
#'
#' Parameters that greatly affect the quality of the layout are the number of nearest
#' neighbors used in the graph (\code{num.nn}) and the aggresiveness of the
#' \code{cut.unconnected.segments} parameter. Additionally, the layout is largely
#' reproducible if given the same starting conditions, but minor changes can affect
#' it quite a bit, so it is sometimes worth varying \code{num.nn} by 1 across a small range
#' (i.e. try 5 layouts with \code{num.nn=120-125}) and choosing the best one. Finally,
#' sometimes different branches can overlap after they have separated, and for a final
#' publication-ready layout, it can make sense to hand-tune the layout some. This can
#' be done using the \code{\link{treeForceRotateCoords}} and \code{\link{treeForceTranslateCoords}}
#' functions to move or rotate sections of the tree. (For instance, in the zebrafish
#' layout, we used these functions to move the completely disconnected EVL and PGC cells
#' into place and to increase the angular distance at the first branchpoint in the blastoderm
#' so that the more fine-grained branches from the different germ layers didn't overlap.)
#'
#' Thanks for Dorde Relic for debugging and creating the \code{remove.duplicate.cells} parameter and Simon Cai for debugging that led to the \code{cell.minimum.walks} parameter.
#'
#' @importFrom RANN nn2
#' @importFrom stats quantile
#' @importFrom reshape2 melt
#' @importFrom igraph graph_from_data_frame layout_with_fr layout_with_drl layout_with_kk vcount V distances is_connected decompose
#'
#' @param object An URD object
#' @param num.nn (Numeric) Number of nearest neighbors to use. (\code{NULL} will use the square root of the number of cells as default.)
#' @param method (Character: "fr", "drl", or "kk") Which force-directed layout algorithm to use (\link[igraph:layout_with_fr]{Fruchterman-Reingold}, \link[igraph:layout_with_drl]{DrL}, or \link[igraph:layout_with_kk]{Kamada-Kawai})
#' @param cells.to.do (Character vector) Cells to use in the layout (default \code{NULL} is all cells in the tree.)
#' @param cell.minimum.walks (Numeric) Minimum number of times a cell must have been visited by random walks in order to be included in the force-directed layout (Cells that have been visited only a few times are more likely to be assigned incorrectly or poorly.)
#' @param cut.outlier.cells (Numeric) If desired, omit cells with unusual second nearest neighbor distances (i.e. those that are likely outliers only well connected to one other cell). Parameter is given as a factor of the interquartile range calculated across all cells' distance to their second nearest neighbor. (Default is not to omit any cells.)
#' @param cut.outlier.edges (Numeric) If desired, cut edges in the nearest neighbor graph with unusually long distances. Parameters is given as a factor of the interquartile range calculated across all edges in the graph. (Default is not to cut any edges based on their length.)
#' @param max.pseudotime.diff (Numeric) If desired, cut edges in the nearest neighbor graph between cells with longer difference in pseudotime. (Default is not to cut any edges based on their pseudotime.)
#' @param cut.unconnected.segments (Numeric) Cut connections in the nearest-neighbor graph to cells that are
#' more segments away in the dendrogram structure. For instance, the most aggressive setting (1) will only
#' maintain that are at most 1 segment away (so it will maintain connections within a segment and to that
#' segment's parent or children. Higher values permit more distant connections. The default value is 2,
#' which permits connections up to two segments away (i.e. connections within a segment, to that segment's
#' parent, grandparent, children, grandchildren, and siblings.) \code{NA} or \code{NULL} disable this setting
#' and permit all connections.
#' @param min.final.neighbors (Numeric) After trimming outlier and unconnected connections in the nearest
#' neighbor graph, remove any cells that remain connected to fewer than this many other cells.
#' @param remove.duplicate.cells (Logical) If cells have exactly the same visitation during the random walks from each tip, this can create problems in the force-directed layout where cells are pushed to the outside of the layout. To avoid this if \code{remove.duplicate.cells=T}, only one cell from each group of cells with duplicated coordinates will be used in the layout.
#' @param tips (Character vector) Tips for which walk visitation data should be used in the construction of the nearest neighbor graph. (Default is all tips )
#' @param coords (Matrix: Cells as rows, 2 columns) Starting coordinates for the force directed layout. Default (\code{"auto"}) takes them from the cell layout of the dendrogram.
#' @param start.temp (Numeric) Starting temperature for the force-directed layout (if \code{method="fr"}), which controls how much cells can move in the initial iterations of the algorithm. Default (\code{NULL}) is the square root of the number of cells.
#' @param n.iter (Numeric or \code{NULL}) Sets the number (or maximum number) of iterations used in the force-directed layout; if \code{NULL}, observes defaults suggested by igraph.
#' @param density.neighbors (Numeric) Distance to this nearest neighbor (default is 10th nearest neighbor) is used as a proxy for local density in the force-directed layout. This is used by \code{\link{plotTreeForce}} if \code{density.alpha=T} for increasing transparency in more high density regions of the layout. This can be adjusted after generating the layout by re-running the \code{\link{fdlDensity}} function.
#' @param plot.outlier.cuts (Logical) If \code{cut.outlier.cells=T} or \code{cut.outlier.edges=T}, this displays the
#' @param verbose (Logical) Print progress and time stamps?
#'
#' @examples
#' object.built <- treeForceDirectedLayout(object.built, num.nn=120, pseudotime="pseudotime", method = "fr", dim = 2, cells.to.do = robustly.visited.cells, tips=final.tips, cut.unconnected.segments = 2, min.final.neighbors=4, verbose=T)
#'
#' @export
treeForceDirectedLayout <- function(object, num.nn=NULL, method=c("fr", "drl", "kk"), cells.to.do=NULL, cell.minimum.walks=1, cut.outlier.cells=NULL, cut.outlier.edges=NULL, max.pseudotime.diff=NULL, cut.unconnected.segments=2, min.final.neighbors=2, remove.duplicate.cells=T, tips=object@tree$tips, coords="auto", start.temp=NULL, n.iter=NULL, density.neighbors=10, plot.outlier.cuts=F, verbose=F) {
# Params
if (length(method) > 1) method <- method[1]
if (is.null(cells.to.do)) cells.to.do <- rownames(object@diff.data)
starting.cells <- length(cells.to.do)
if (is.null(num.nn)) num.nn <- ceiling(sqrt(length(cells.to.do)))
if (is.null(tips)) tips <- object@tree$tips
if (class(coords)=="character" && coords=="auto") {
coords <- as.matrix(object@tree$cell.layout[cells.to.do,c("x","y")])
} else if (!is.null(coords) && class(coords) != "matrix") {
stop("coords must either be 'auto', NULL, or a matrix.")
}
dim=2
if (verbose) print(paste(Sys.time(), ": Starting with parameters", method, num.nn, "NN", dim, "D", length(cells.to.do), "cells"))
# Get rid of any cells that don't have a pseudotime (or, if cut.unconnected.segments is turned on, any cells that aren't part of the tree)
cells.no.pseudotime <- cells.to.do[which(is.na(object@tree$pseudotime[cells.to.do]))]
if (!is.na(cut.unconnected.segments) & !is.null(cut.unconnected.segments) & cut.unconnected.segments > 0) {
cells.not.in.tree <- setdiff(cells.to.do, unlist(object@tree$cells.in.segment))
if (verbose) print(paste0("Removing ", length(unique(c(cells.no.pseudotime, cells.not.in.tree))), " cells that are not assigned a pseudotime or a segment in the tree."))
cells.to.do <- setdiff(cells.to.do, c(cells.no.pseudotime, cells.not.in.tree))
} else {
if (verbose) print(paste0("Removing ", length(cells.no.pseudotime), " cells that are not assigned a pseudotime."))
cells.to.do <- setdiff(cells.to.do, cells.no.pseudotime)
}
# Get and normalize walk data, then add pseudotime (unless NULL in which case don't use)
if (verbose) print(paste0(Sys.time(), ": Preparing walk data."))
walk.data <- object@diff.data[cells.to.do,paste0("visitfreq.raw.", object@tree$tips)]
walk.data$pseudotime <- object@tree$pseudotime[cells.to.do]
walk.total <- apply(walk.data, 1, sum)
# If some cells weren't walked enough
if (any(walk.total < cell.minimum.walks)) {
if (verbose) print(paste0("Removing ", length(which(walk.total < cell.minimum.walks)), " cells that were visited fewer than ", cell.minimum.walks, " times by random walks."))
# Remove those cells
walk.data <- walk.data[which(walk.total >= cell.minimum.walks),]
# And recalculate walk totals
walk.total <- apply(walk.data, 1, sum)
}
walk.data <- sweep(walk.data, 1, walk.total, "/")
walk.data <- as.matrix(walk.data)
duped <- which(duplicated(walk.data))
if (length(duped) > 0 && remove.duplicate.cells) {
warning(length(duped), " cells have duplicate random walk coordinates and are being removed from the layout.")
walk.data <- unique(walk.data)
} else if (length(duped) > 0) {
warning(length(duped), " cells have duplicate random walk coordinates. This may cause a problem in the layout. If so, set remove.duplicate.cells=T")
}
# Calculate nearest neighbor graph
if (verbose) print(paste0(Sys.time(), ": Calculating nearest neighbor graph."))
walk.nn <- RANN::nn2(data=walk.data,query=walk.data,k=max(num.nn)+1, treetype = "kd", searchtype="priority")
rownames(walk.nn$nn.idx) <- rownames(walk.data)
walk.nn$nn.idx <- walk.nn$nn.idx[,2:(num.nn+1)]
rownames(walk.nn$nn.dists) <- rownames(walk.data)
walk.nn$nn.dists <- walk.nn$nn.dists[,2:(num.nn+1)]
walk.nn$nn.label <- apply(walk.nn$nn.idx, 2, function(y) rownames(walk.data)[y])
rownames(walk.nn$nn.label) <- rownames(walk.data)
# Cull outlier cells with unusual second nearest neighbor distances
if (!is.null(cut.outlier.cells)) {
q <- quantile(walk.nn$nn.dists[,2])
outer.fence <- q[4] + cut.outlier.cells*(q[4]-q[1])
if (plot.outlier.cuts) {
hist(walk.nn$nn.dists[,2], breaks=100, main="Outlier Removal: Second Nearest Neighbor", xlab="Distance to 2nd Nearest Neighbor")
abline(v=outer.fence, col='red')
}
if (verbose) print(paste0(Sys.time(), ": Removing ", round(length(which(walk.nn$nn.dists[,2] > outer.fence)) / starting.cells, digits=2), "% of cells as outliers."))
cells.keep <- which(walk.nn$nn.dists[,2] <= outer.fence)
} else {
cells.keep <- 1:dim(walk.data)[1]
}
# Turn nearest neighbors into edge list
if (verbose) print(paste0(Sys.time(), ": Preparing edge list."))
edges <- reshape2::melt(walk.nn$nn.label[cells.keep,], stringsAsFactors=F)
dists <- reshape2::melt(walk.nn$nn.dists[cells.keep,])
edges <- edges[,c(1,3)]
names(edges) <- c("V1", "V2")
edges$dists <- dists$value
edges$V1 <- as.character(edges$V1)
edges$V2 <- as.character(edges$V2)
starting.edges <- dim(edges)[1]
# Cut connections with unusually long distances
if (!is.null(cut.outlier.edges)) {
q <- quantile(edges$dists)
outer.fence <- q[4] + cut.outlier.edges*(q[4]-q[1])
if (verbose) print(paste0(Sys.time(), ": Removing ", round(length(which(edges$dists > outer.fence)) / starting.edges * 100, digits=3), "% of edges as outliers."))
if (plot.outlier.cuts) {
hist(edges$dists, breaks=100, main="Outlier Removal: Edge distances", xlab="Edge distance")
abline(v=outer.fence, col='red')
}
edges <- edges[edges$dists <= outer.fence,]
}
# Cut connections with overly large pseudotime differences
if (!is.null(max.pseudotime.diff)) {
edges$pt1 <- object@tree$pseudotime[edges$V1]
edges$pt2 <- object@tree$pseudotime[edges$V2]
edges$dpt <- abs(edges$pt1 - edges$pt2)
if (verbose) print(paste0(Sys.time(), ": Removing ", round(length(which(edges$dpt > max.pseudotime.diff)) / starting.edges * 100, digits=3), "% of edges with too large pseudotime difference."))
edges <- edges[which(edges$dpt <= max.pseudotime.diff),]
}
# Cut edges between overly distant (e.g. non-connected) segments
if (!is.na(cut.unconnected.segments) & !is.null(cut.unconnected.segments) & cut.unconnected.segments > 0) {
# Determine which segments each connected cell belongs to
edges$seg1 <- object@diff.data[edges$V1, "segment"]
edges$seg2 <- object@diff.data[edges$V2, "segment"]
# Determine distance in the tree between each segment
seg.dist <- igraph::distances(object@tree$tree.igraph, mode="all")
edges$seg.dist <- apply(edges[,c("seg1","seg2")], 1, function(segs) seg.dist[segs[1],segs[2]])
if (verbose) print(paste0(Sys.time(), ": Removing ", round(length(which(edges$seg.dist > cut.unconnected.segments)) / starting.edges * 100, digits=2), "% of edges that are between segments with distance > ", cut.unconnected.segments))
edges <- edges[edges$seg.dist <= cut.unconnected.segments,]
}
# Trim cells that are no longer well connected to the graph
if (verbose) print(paste0(Sys.time(), ": Trimming cells that are no longer well connected."))
connections.remaining <- table(edges[edges$dists > 0,"V1"])
cells.with.enough.connections <- names(which(connections.remaining >= min.final.neighbors))
cells.without.enough.connections <- length(unique(edges$V1)) - length(cells.with.enough.connections)
while(cells.without.enough.connections > 0) {
edges <- edges[which(edges$V1 %in% cells.with.enough.connections & edges$V2 %in% cells.with.enough.connections),]
connections.remaining <- table(edges[edges$dists > 0,"V1"])
cells.with.enough.connections <- names(which(connections.remaining >= min.final.neighbors))
cells.without.enough.connections <- length(unique(edges$V1)) - length(cells.with.enough.connections)
}
if (verbose) print(paste0(Sys.time(), ": ", round(length(cells.with.enough.connections)/starting.cells*100, digits=2), "% of starting cells preserved."))
# Create igraph representation
if (verbose) print(paste0(Sys.time(), ": Preparing igraph object."))
if (method=="kk") {
# kk requires that weights be higher for points that should be far apart -- use distance directly
edges$weight <- edges$dists
} else {
# everything else places points with high weight nearby -- use inverse of distance as graph weight
edges$weight <- 1/edges$dists
}
object@tree$walks.force.edges <- edges
igraph.walk.weights <- igraph::graph_from_data_frame(edges, directed=F)
if (!is.null(coords)) {
cells.remain <- unique(c(edges$V1, edges$V2))
coords.remain <- rownames(coords) %in% cells.remain
coords <- coords[coords.remain,]
}
# DrL force-directed layout only works on connected graphs, so check for graph connectivity
if (method == "drl" && !igraph::is_connected(igraph.walk.weights)) {
igraph.walk.weights.decomposed <- decompose(igraph.walk.weights)
graph.keep <- which.max(unlist(lapply(igraph.walk.weights.decomposed, function(x) length(V(x)))))
warning("DrL method can only operate on fully connected graphs. Using more nearest neighbors will increase graph connectivity. For now, using largest subgraph: ", length(V(igraph.walk.weights.decomposed[[graph.keep]])), " of ", length(V(igraph.walk.weights)), " cells.")
igraph.walk.weights <- igraph.walk.weights.decomposed[[graph.keep]]
coords <- coords[rownames(coords) %in% names(V(igraph.walk.weights)),]
cells.with.enough.connections <- intersect(cells.with.enough.connections, names(V(igraph.walk.weights)))
}
# Define start temp, if set to default
if (is.null(start.temp)) start.temp <- sqrt(igraph::vcount(igraph.walk.weights))
# Define number of iterations to do
if (is.null(n.iter)) {
if (method == "fr") n.iter = 500
if (method == "kk") n.iter = 50*igraph::vcount(igraph.walk.weights)
}
# Do force-directed layout
if (verbose) print(paste0(Sys.time(), ": Doing force-directed layout."))
if (method=="fr") igraph.walk.layout <- igraph::layout_with_fr(igraph.walk.weights, dim=dim, coords=coords, start.temp=start.temp, niter=n.iter)
if (method=="drl") igraph.walk.layout <- igraph::layout_with_drl(igraph.walk.weights, dim=dim, options=list(edge.cut=0))
if (method=="kk") igraph.walk.layout <- igraph::layout_with_kk(igraph.walk.weights, dim=dim, coords=coords, maxiter=n.iter)
# Store the layout
if (dim==2) {
# Calculate weighted pseudotime of cells' neighbors for tree layout.
if (verbose) print(paste0(Sys.time(), ": Calculating Z."))
neighbor.pt <- unlist(lapply(cells.with.enough.connections, function(i) {
# Deal with points with 0 distance, by setting their weight to the sum of all other weights.
weights <- 1/walk.nn$nn.dists[i,]
weights[!is.finite(weights)] <- sum(weights[is.finite(weights)])
weighted.mean(x=object@tree$pseudotime[walk.nn$nn.label[i,]], w=weights)
}))
names(neighbor.pt) <- cells.with.enough.connections
# Store the layout
object@tree$walks.force.layout <- as.data.frame(igraph.walk.layout, stringsAsFactors=F)
names(object@tree$walks.force.layout) <- c("x","y")
rownames(object@tree$walks.force.layout) <- igraph::V(igraph.walk.weights)$name
# Remove points with NA values for force-directed layout
fdl.failed <- length(which(!(complete.cases(object@tree$walks.force.layout))))
if (fdl.failed > 0) {
warning(fdl.failed, " cells were not assigned coordinates during the force-directed layout.")
object@tree$walks.force.layout <- object@tree$walks.force.layout[complete.cases(object@tree$walks.force.layout),]
}
# Normalize neighbor.pt range
neighbor.pt.factor <- mean(apply(object@tree$walks.force.layout[,c("x","y")], 2, max)) / max(neighbor.pt, na.rm=T)
object@tree$walks.force.layout$telescope.pt <- neighbor.pt[rownames(object@tree$walks.force.layout)] * neighbor.pt.factor
# Calculate local density for adjusting alpha
if (verbose) print(paste0(Sys.time(), ": Calculating local density."))
object <- fdlDensity(object, neighbor=density.neighbors)
} else if (dim==3) {
object@tree$walks.force.layout.3d <- as.data.frame(igraph.walk.layout, stringsAsFactors=F)
names(object@tree$walks.force.layout.3d) <- c("x","y","z")
rownames(object@tree$walks.force.layout.3d) <- igraph::V(igraph.walk.weights)$name
}
# Position labels if segments have been named.
if (!is.null(object@tree$segment.names)) {
object@tree$walks.force.labels <- treeForcePositionLabels(object)
}
# Return the object
if (verbose) print(paste0(Sys.time(), ": Finished."))
return(object)
}
#' Plot tree force-directed layout in 2D
#'
#' @import ggplot2
#'
#' @param object An URD object
#'
#' @export
plotTreeForce2D <- function(object, label=NULL, label.type="search", title=label, show.points=T, point.alpha=1, point.size=1, show.neighbors=F, neighbors.max=10000, colors=NULL) {
# Load colors if not specified
if (is.null(colors)) colors <- defaultURDContinuousColors()
gg.data <- object@tree$walks.force.layout
if (!is.null(label)) color.data <- data.for.plot(object = object, label = label, label.type = label.type, as.color = F, as.discrete.list=T, cells.use = rownames(gg.data))
gg.data$expression <- color.data$data
gg <- ggplot() + theme_bw() + labs(x="", y="", color="", title=title)
if (show.neighbors) {
edge.data <- object@tree$walks.force.edges
if (dim(edge.data)[1] > neighbors.max) edge.data <- edge.data[sample(1:dim(edge.data)[1], neighbors.max, replace=F),]
edge.data[,c("x1","y1")] <- gg.data[edge.data$V1,c("x","y")]
edge.data[,c("x2","y2")] <- gg.data[edge.data$V2,c("x","y")]
if (show.points) gg <- gg + geom_segment(data=edge.data, aes(x=x1, y=y1, xend=x2, yend=y2), color='black', alpha=0.02)
if (!show.points) gg <- gg + geom_segment(data=edge.data, aes(x=x1, y=y1, xend=x2, yend=y2, color=weight), alpha=0.04)
}
if (show.points) gg <- gg + geom_point(data=gg.data, aes(x=x, y=y, color=expression), size=point.size, alpha=point.alpha)
if (!color.data$discrete) gg <- gg + scale_color_gradientn(colors = colors)
return(gg)
}
#' Plot force-directed layout of tree
#'
#' This plots data on the force-directed layout of an URD tree, generated by
#' \code{\link{treeForceDirectedLayout}}. Once the tree is rotated into a
#' perspective that is appealing, \code{\link{plotTreeForceStore3DView}} can be
#' used to store that view in the URD object for additional plots.
#'
#' The transparency of points is typically adjusted in two ways: (1) Points with
#' low expression are made more transparent, which is controlled through the
#' \code{fade.below} and \code{alpha.fade} parameters. (2) Points in higher
#' density regions of the plot are made more transparent, which is controlled by
#' the \code{density.alpha} parameter (and the local density is determined by
#' the \code{density.neighbors} parameter of \code{\link{treeForceDirectedLayout}},
#' and can be adjusted after the fact by re-running \code{\link{fdlDensity}}).
#'
#' Plots can be easily saved using \code{\link{rgl:rgl.snapshot}} and plots
#' can be closed with \code{\link{rgl:rg.close}} if plotting many genes
#' using a loop.
#'
#' @param object An URD object
#' @param label (Character) Data to use for coloring tree (see \link{data.for.plot})
#' @param label.type (Character) Where to find data for label (default \code{"search"} auto-detects, see \link{data.for.plot})
#' @param view (Character) Name of view to use (i.e. the name of an entry in \code{@@tree$force.view.list}) created by \link{plotTreeForceStore3DView}. Default is last view created.
#' @param alpha (Numeric) Maximum transparency of points
#' @param alpha.fade (Numeric) Minimum transparency of points
#' @param size (Numeric) Size of points in the plot
#' @param title (Character) Title to add to the plot. (This is sensitive to resizing the window after plotting, but if a view is stored, the window will be resized before the title is added, so it will be acceptable resolution for figures.)
#' @param title.cex (Numeric) Adjust the title font size
#' @param title.line (Numeric) Adjust the position of the title. Positive numbers move the title upward.
#' @param label.tips (Logical) Should text identifying the tips of the tree be placed in the force directed layout? Defaults to TRUE if \code{@@tree$segment.names} has been set.
#' @param use.short.names (Logical) Should short names from \code{@@tree$segment.names.short} be used? Defaults to TRUE if those values have been set.
#' @param seg.show (Character vector) Segments of the tree to put on the plot (Default \code{NULL} is all segments)
#' @param cells.show (Character vector) Cells of the tree to show (Default \code{NULL} is all cells)
#' @param fade.below (Numeric) If desired, transparency can be lowered for points with low expression values (thereby highlighting the portions of the tree with expression). This value is the portion of the range of expression values to fade. (The default, 2/9, fades the bottom two-ninths of the expression values.) A fully 'faded' point will have alpha of \code{alpha.fade}.
#' @param density.alpha (Logical) Should points with short distance to nearest neighbors (i.e. in denser regions of the plot) be more transparent?
#' @param label.spacing (Numeric) How far should labels be from the detected tips
#' @param text.cex (Numeric) Size of the label text
#' @param colors (Character vector) Vector of colors to use if plotting continuous data
#' @param discrete.colors (Character vector) Vector of colors to use if plotting discrete data
#' @param color.limits (Numeric vector, length 2) Minimum and maximum values for color scale. Default \code{NULL} auto-detects.
#' @param symmetric.color.scale (Logical) Should the color scale be symmetric and centered around 0? (Default \code{NULL} is \code{FALSE} if all values are positive, and \code{TRUE} if both positive and negative values are present.)
#'
#' @return Nothing. Produces a plot using the \code{rgl} package, displayed in an X11 window.
#'
#' @examples
#' # Define the colors used below
#' stage.colors <- c("#CCCCCC", RColorBrewer::brewer.pal(9, "Set1")[9], RColorBrewer::brewer.pal(12, "Paired")[c(9,10,7,8,5,6,3,4,1,2)])
#' fire.with.grey <- c("#CECECE", "#DDC998", RColorBrewer::brewer.pal(9, "YlOrRd")[3:9])
#' branch.colors <- c("#CECECE", "#E6298B")
#'
#' # Plot a discrete characteristic like stage or clustering results
#' plotTreeForce(object, "stage.nice", title="STAGE", title.line=1, discrete.colors = stage.colors, alpha=0.4)
#'
#' # Plot gene expression
#' plotTreeForce(object, "NOTO", title="NOTO expression", title.line=1, colors = fire.with.grey)
#'
#' Define a 'clustering' based on whether or not cells are in a given trajectory then plot it
#' object <- groupFromCells(object, group.id="lineage_Tailbud", cells=cellsAlongLineage(object, "Tailbud", remove.root=F))
#' plotTreeForce(object, "lineage_Tailbud", title="Tailbud", title.line=1, discrete.colors=branch.colors, alpha=0.4)
#'
#' @export
plotTreeForce <- function(object, label, label.type="search", view="default", alpha=0.8, alpha.fade=0.1, size=5, title=NULL, title.cex=3, title.line=0, label.tips=(!is.null(object@tree$segment.names) | !is.null(object@tree$segment.names.short)), use.short.names=!is.null(object@tree$segment.names.short), seg.show=NULL, cells.show=NULL, fade.below=(2/9), density.alpha=T, label.spacing=5, text.cex=0.8, colors=NULL, discrete.colors=NULL, color.limits=NULL, symmetric.color.scale=NULL) {
if (requireNamespace("rgl", quietly = TRUE)) {
# Figure out which view to use
if (!is.null(view)) {
if (view=="default") {
if ("force.view.default" %in% names(object@tree)) view <- object@tree$force.view.default
else view <- NULL
}
view <- object@tree$force.view.list[[view]]
}
# Get layout data
gg.data <- object@tree$walks.force.layout
gg.data$segment <- object@diff.data[rownames(gg.data),"segment"]
# Get expression data, use to determine parameters of color scale, and get default colors
expression.data <- data.for.plot(object=object, label=label, label.type=label.type, as.color=F, cells.use=rownames(gg.data), as.discrete.list=T)
if (!expression.data$discrete && is.null(symmetric.color.scale)) symmetric.color.scale <- min(expression.data$data) < 0
if (!expression.data$discrete && is.null(color.limits)) {
if (symmetric.color.scale) {
color.sv <- max(abs(expression.data$data))
color.limits <- c(-1*color.sv, color.sv)
} else {
color.limits <- c(0, max(expression.data$data))
}
}
if (is.null(colors) && !expression.data$discrete) colors <- defaultURDContinuousColors(symmetric = symmetric.color.scale)
# Set color, alpha, and point size
if (!is.null(label)) {
color.data <- data.for.plot(object = object, label = label, label.type = label.type, as.color = T, as.discrete.list=T, cells.use = rownames(gg.data), continuous.colors=colors, colors.use = discrete.colors, continuous.color.limits = color.limits)
gg.data$expression <- color.data$data
}
gg.data$alpha <- alpha
gg.data$size <- size
# Focus on expression
# For this, grey out any thing below focus, set alpha low for those, and increase progressively until 2*focus.
if (fade.below > 0 & !color.data$discrete) {
# Figure out ranges
expression.range <- range(abs(expression.data$data))
fade.range <- diff(expression.range) * fade.below + expression.range[1]
# Adjust gg.data
fade <- (fade.range-abs(expression.data$data)+expression.range[1])/fade.range
fade[fade < 0] <- 0
gg.data$alpha <- alpha - ((alpha-alpha.fade) * fade)
}
# Open 3D view
if (!is.null(view)) {
rgl::open3d(
zoom=view$rgl.setting$zoom,
scale=view$rgl.setting$scale,
userMatrix=view$rgl.setting$userMatrix,
windowRect=view$rgl.setting$windowRect
)
} else {
rgl::open3d()
}
# Adjust transparency for local density
if (density.alpha) {
dr <- range(object@tree$walks.force.layout$n.dist)
n.dist <- object@tree$walks.force.layout$n.dist / dr[2]
density.adj <- sqrt(n.dist)
density.adj <- density.adj / median(density.adj)
gg.data$alpha <- gg.data$alpha * density.adj
}
# Convert alpha to 0 on cells not to show (so the plot size doesn't get changed)
if (!is.null(seg.show) & is.null(cells.show)) {
segs.show <- segChildrenAll(object, seg.show, include.self=T)
cells.show <- rownames(gg.data)[which(gg.data$segment %in% segs.show)]
}
if (!is.null(cells.show)) {
cells.hide <- setdiff(rownames(gg.data), cells.show)
gg.data[cells.hide,"alpha"] <- 0
}
# Add the points to the plot
rgl::points3d(x=gg.data$x, y=gg.data$y, z=gg.data$telescope.pt, col=gg.data$expression, alpha=gg.data$alpha, size=size)
# Add title to plot
if (!is.null(title)) {
Sys.sleep(0.1)
rgl::bgplot3d({plot.new(); title(main=title, line=title.line, cex.main=title.cex)})
Sys.sleep(0.1)
}
# Add labels to tips
if (label.tips) {
# Check for existing tree labels
if (!is.null(object@tree$walks.force.labels)) {
tip.labels <- object@tree$walks.force.labels
} else {
# If there aren't any, generate some.
tip.labels <- treeForcePositionLabels(object)
}
# Choose short or long labels
if (use.short.names) tip.labels$label <- tip.labels$name.short else tip.labels$label <- tip.labels$name
# Limit to labels of cells that appear in the plot
segs.to.label <- intersect(unique(gg.data[which(gg.data$alpha > 0),"segment"]), tip.labels$seg)
tip.labels <- tip.labels[which(tip.labels$seg %in% segs.to.label),]
# Add the text
rgl::text3d(x=tip.labels$x, y=tip.labels$y, z=tip.labels$z, text=tip.labels$label, adj=0.5, cex=text.cex)
}
# Add a brief pause to ensure that rendering completes before moving on to next function.
Sys.sleep(0.1)
} else {
stop("Package rgl is required for this function. To install: install.packages('rgl')\n")
}
}
#' Plot force-directed layout of tree
#'
#' This plots data on the force-directed layout of an URD tree, generated by
#' \code{\link{treeForceDirectedLayout}}, using red and green (which become yellow
#' when they overlap to display two pieces of data simultaneously. See
#' \code{\link{plotTreeForce}} for more details.
#'
#' @param object An URD object
#' @param label.red (Character) Data to use for coloring tree (red channel) (see \link{data.for.plot})
#' @param label.green (Character) Data to use for coloring tree (green channel) (see \link{data.for.plot})
#' @param label.type.red (Character) Where to find data for red channel label (default \code{"search"} auto-detects, see \link{data.for.plot})
#' @param label.type.green (Character) Where to find data for green channel label (default \code{"search"} auto-detects, see \link{data.for.plot})
#' @param view (Character) Name of view to use (i.e. the name of an entry in \code{@@tree$force.view.list}) created by \link{plotTreeForceStore3DView}. Default is last view created.
#' @param alpha (Numeric) Maximum transparency of points
#' @param alpha.fade (Numeric) Minimum transparency of points
#' @param size (Numeric) Size of points in the plot
#' @param title (Character) Title to add to the plot. (This is sensitive to resizing the window after plotting, but if a view is stored, the window will be resized before the title is added, so it will be acceptable resolution for figures.)
#' @param title.cex (Numeric) Adjust the title font size
#' @param title.line (Numeric) Adjust the position of the title. Positive numbers move the title upward.
#' @param label.tips (Logical) Should text identifying the tips of the tree be placed in the force directed layout? Defaults to TRUE if \code{@@tree$segment.names} has been set.
#' @param use.short.names (Logical) Should short names from \code{@@tree$segment.names.short} be used? Defaults to TRUE if those values have been set.
#' @param seg.show (Character vector) Segments of the tree to put on the plot (Default \code{NULL} is all segments)
#' @param cells.show (Character vector) Cells of the tree to show (Default \code{NULL} is all cells)
#' @param fade.below (Numeric) If desired, transparency can be lowered for points with low expression values (thereby highlighting the portions of the tree with expression). This value is the portion of the range of expression values to fade. (The default, 1/3, fades the bottom third of the expression values.) A fully 'faded' point will have alpha of \code{alpha.fade}.
#' @param fade.method (Character: "max", "min", or "mean") Should fade be calculated on max expression of the two labels (emphasize any cells expressing either label), min expression of the two labels (emphasize only cells co-expressing both labels), or mean expression
#' @param density.alpha (Logical) Should points with short distance to nearest neighbors (i.e. in denser regions of the plot) be more transparent?
#' @param label.spacing (Numeric) How far should labels be from the detected tips
#' @param text.cex (Numeric) Size of the label text
#'
#' @return Nothing. Produces a plot using the \code{rgl} package, displayed in an X11 window.
#'
#' @examples
#' # Focus on all cells that express either gene
#' plotTreeForceDual(object, label.red="SOX17", label.green="SOX32", title="SOX32 (green) + SOX17 (red)")
#'
#' # Focus on only the co-expressing cells by using fade.method="min"
#' plotTreeForceDual(object, label.red="SOX17", label.green="SOX32", title="SOX32 (green) + SOX17 (red)", fade.method="min")
#'
#' @export
plotTreeForceDual <- function(object, label.red, label.green, label.red.type="search", label.green.type="search", view="default", alpha=0.8, alpha.fade=0.025, size=5, title=NULL, title.cex=3, title.line=0, label.tips=(!is.null(object@tree$segment.names) | !is.null(object@tree$segment.names.short)), use.short.names=!is.null(object@tree$segment.names.short), seg.show=NULL, cells.show=NULL, fade.below=(1/3), fade.method=c("max", "min", "mean"), density.alpha=T, label.spacing=5, text.cex=0.8) {
if (requireNamespace("rgl", quietly = TRUE)) {
# Figure out which view to use
if (!is.null(view)) {
if (view=="default") {
if ("force.view.default" %in% names(object@tree)) view <- object@tree$force.view.default
else view <- NULL
}
view <- object@tree$force.view.list[[view]]
}
# Get layout data
gg.data <- object@tree$walks.force.layout
gg.data$segment <- object@diff.data[rownames(gg.data),"segment"]
# Get expression data
red.data <- data.for.plot(object=object, label=label.red, label.type=label.red.type, as.single.color = T, as.discrete.list=T, cells.use=rownames(gg.data))
green.data <- data.for.plot(object=object, label=label.green, label.type=label.green.type, as.single.color = T, as.discrete.list=T, cells.use=rownames(gg.data))
if (red.data$discrete | green.data$discrete) stop("For dual-color plots, both labels must be continuous values.")
gg.data$expression <- rgb(red.data$data, green.data$data, 0)
gg.data$alpha <- alpha
gg.data$size <- size
# Focus on expression
# For this, grey out any thing below focus, set alpha low for those, and increase progressively until 2*focus.
if (fade.below > 0) {
# Figure out value for fading
if (length(fade.method) > 1) fade.method <- fade.method[1]
if (tolower(fade.method) == "max") {
# Use max of red or green for fading
expression.data <- pmax(red.data$data, green.data$data)
} else if (tolower(fade.method) == "min") {
# Use min of red or green for fading
expression.data <- pmin(red.data$data, green.data$data)
} else if (tolower(fade.method) == "mean") {
# Use min of red or green for fading
expression.data <- apply(data.frame(red.data$data, green.data$data), 1, mean)
} else {
stop("fade.method must be max, min, or mean.")
}
# Adjust gg.data
fade <- (fade.below - expression.data)/fade.below
fade[fade < 0] <- 0
gg.data$alpha <- alpha - ((alpha-alpha.fade) * fade)
}
# Open 3D view
if (!is.null(view)) {
rgl::open3d(
zoom=view$rgl.setting$zoom,
scale=view$rgl.setting$scale,
userMatrix=view$rgl.setting$userMatrix,
windowRect=view$rgl.setting$windowRect
)
} else {
rgl::open3d()
}
# Adjust transparency for local density
if (density.alpha) {
dr <- range(object@tree$walks.force.layout$n.dist)
n.dist <- object@tree$walks.force.layout$n.dist / dr[2]
density.adj <- sqrt(n.dist)
density.adj <- density.adj / median(density.adj)
gg.data$alpha <- gg.data$alpha * density.adj
}
# Convert alpha to 0 on cells not to show (so the plot size doesn't get changed)
if (!is.null(seg.show) & is.null(cells.show)) {
segs.show <- segChildrenAll(object, seg.show, include.self=T)
cells.show <- rownames(gg.data)[which(gg.data$segment %in% segs.show)]
}
if (!is.null(cells.show)) {
cells.hide <- setdiff(rownames(gg.data), cells.show)
gg.data[cells.hide,"alpha"] <- 0
}
# Add the points to the plot
rgl::points3d(x=gg.data$x, y=gg.data$y, z=gg.data$telescope.pt, col=gg.data$expression, alpha=gg.data$alpha, size=size)
# Add title to plot
if (!is.null(title)) {
Sys.sleep(0.1)
rgl::bgplot3d({plot.new(); title(main=title, line=title.line, cex.main=title.cex)})
Sys.sleep(0.1)
}
# Add labels to tips
if (label.tips) {
# Check for existing tree labels
if (!is.null(object@tree$walks.force.labels)) {
tip.labels <- object@tree$walks.force.labels
} else {
# If there aren't any, generate some.
tip.labels <- treeForcePositionLabels(object)
}
# Choose short or long labels
if (use.short.names) tip.labels$label <- tip.labels$name.short else tip.labels$label <- tip.labels$name
# Limit to labels of cells that appear in the plot
segs.to.label <- intersect(unique(gg.data[which(gg.data$alpha > 0),"segment"]), tip.labels$seg)
tip.labels <- tip.labels[which(tip.labels$seg %in% segs.to.label),]
# Add the text
rgl::text3d(x=tip.labels$x, y=tip.labels$y, z=tip.labels$z, text=tip.labels$label, adj=0.5, cex=text.cex)
}
# Add a brief pause to ensure that rendering completes before moving on to next function.
Sys.sleep(0.1)
} else {
stop("Package rgl is required for this function. To install: install.packages('rgl')\n")
}
}
#' Determine position for segment labels in force-directed layout
#'
#' Finds cells near the end of each tip, and attempts to calculate a vector from
#' them to determine a good location for a label near each tip. Called automatically
#' by \code{\link{treeForceDirectedLayout}}, but can be re-run if needed.
#'
#' @param object An URD object
#' @param label.spacing (Numeric) Length of vector from final cell in each tip to label
#'
#' @return A data.frame with rows as tips, and columns containing labels and coordinates.
#' This is normally stored in \code{@@tree$walks.force.labels}.
#'
#' @examples
#' object@tree$walks.force.labels <- treeForcePositionLabels(object, label.spacing=2) # Move labels closer
#'
#' @export
treeForcePositionLabels <- function(object, label.spacing=5) {
# Grab the tips and their names
tip.labels <- data.frame(
seg=segTerminal(object),
name=object@tree$segment.names[segTerminal(object)],
name.short=object@tree$segment.names.short[segTerminal(object)],
stringsAsFactors=F, row.names=segTerminal(object)
)
# Anything that wasn't assigned a name previously: NA -> blank
tip.labels[is.na(tip.labels)] <- ""
# Get force-directed layout and add segment and pseudotime information
fdl.data <- object@tree$walks.force.layout
fdl.data$segment <- object@diff.data[rownames(fdl.data),"segment"]
fdl.data$pseudotime <- object@tree$pseudotime[rownames(fdl.data)]
# Order fdl data by pseudotime
fdl.data <- fdl.data[order(fdl.data$pseudotime, decreasing=T),]
# Figure out the location for each label by projecting out the vector at the end of each tip
tip.labels[,c("x","y","z")] <- t(as.data.frame(lapply(tip.labels$seg, function(tip) {
# Get cells from this tip
this.fdl <- fdl.data[which(fdl.data$segment==tip),]
# Find the vector describing direction that the cells at the tip are moving
tip.vec <- as.numeric(apply(this.fdl[c(1,5),c("x","y","telescope.pt")], 2, diff))
# Normalize the length of vector
s <- sqrt(label.spacing^2/(tip.vec[1]^2+tip.vec[2]^2+tip.vec[3]^2))
tip.vec <- s*tip.vec
# Place the label
label.coords <- apply(rbind(this.fdl[1,c("x","y","telescope.pt")], tip.vec),2,sum)
return(label.coords)
})))
return(tip.labels)
}
#' Calculate force-directed layout local density
#'
#' Add distance to nth nearest neighbor to force directed layout so that you can later
#' calculate local density and adjust alpha values.
#'
#' @importFrom RANN nn2
#'
#' @param object An URD object
#' @param neighbor (Numeric) Distance to \code{neighbor}th nearest neighbor is used for
#' density adjustment of point transparency on force-directed layout plots.
fdlDensity <- function(object, neighbor=10, knn=NULL) {
knn <- RANN::nn2(object@tree$walks.force.layout[,c("x","y","telescope.pt")], k = neighbor, treetype = 'kd')
object@tree$walks.force.layout$n.dist <- knn$nn.dists[,neighbor]
return(object)
}
#' Store a 3D view for force-directed layouts
#'
#' The last stored view will be the "default."
#'
#' @param object An URD object
#' @param view.name (Character) Name to use for this view
#'
#' @return An URD object with the parameters for that view stored in \code{@@tree$force.view.list}.
#'
#' @export
plotTreeForceStore3DView <- function(object, view.name) {
if (requireNamespace("rgl", quietly = TRUE)) {
# If not appending to a views list, make one
if (is.null(object@tree$force.view.list)) object@tree$force.view.list <- list()
# Grab rgl plot settings, if desired
rgl.settings <- rgl::par3d(no.readonly=F)
# Make a list of everything
this.view <- list(rgl.settings=rgl.settings)
# Store if in the views list
object@tree$force.view.list[[view.name]] <- this.view
# Store 'default' view
object@tree$force.view.default <- view.name
# Return the final views list
return(object)
} else {
stop("Package rgl is required for this function. To install: install.packages('rgl')\n")
}
}
#' Rotate coordinates in 3D
#'
#' @param m (Matrix) n cells x 3 matrix of coordinates
#' @param a (Numeric) Angle to rotate
#' @param axis (Character) Which axis to rotate around
#'
#' @return n x 3 matrix of rotated coordinates
#'
#' @keywords internal
rotateCoords3d <- function(m, a, axis=c("x","y","z")) {
if (length(axis) > 1) axis <- axis[1]
if (dim(m)[2] != 3) stop("Supply an n x 3 matrix.")
if (axis == "x") r.mat <- matrix(c(1, 0, 0, 0, cos(a), -1*sin(a), 0, sin(a), cos(a)), ncol=3, byrow=T)
if (axis == "y") r.mat <- matrix(c(cos(a), 0, sin(a), 0, 1, 0, -1*sin(a), 0, cos(a)), ncol=3, byrow=T)
if (axis == "z") r.mat <- matrix(c(cos(a), -1*sin(a), 0, sin(a), cos(a), 0, 0, 0, 1), ncol=3, byrow=T)
m %*% r.mat
}
#' Translate coordinates in 3D
#'
#' @param m (Matrix) n cells x 3 coordinates matrix
#' @param x (Numeric) Distance to move in x
#' @param y (Numeric) Distance to move in y
#' @param z (Numeric) Distance to move in z
#'
#' @return n x 3 matrix of translated coordinates
#'
#' @keywords internal
translateCoords3d <- function(m, x, y, z) {
if (dim(m)[2] != 3) stop("Supply an n x 3 matrix.")
sweep(m, 2, c(x, y, z), "+")
}
#' Rotate points in force-directed layout
#'
#' This function can be used in addition to \code{\link{treeForceTranslateCoords}}
#' in order to fine-tune the presentation of a force-directed layout. For instance,
#' this can be used to achieve an improved 2D visualization by modifying overlapping
#' portions of the layout.
#'
#' @param object An URD object
#' @param cells (Character vector) Cells to modify (Default \code{NULL} is all cells in the force-directed layout)
#' @param seg (Character) Instead of specifying cells, just grab all cells from this segment and downstream. Ignored if \code{cells} is specified.
#' @param angle (Numeric) Angle to rotate
#' @param axis (Character: "x", "y", "z") Axis around which to rotate
#' @param around.cell (Character or Numeric) If a character, then the name of a cell to rotate around. If numeric, then will choose the Nth cell by pseudotime. (This can be helpful for rotating a segment around its base to spread segments of the tree apart more. If \code{NULL}, then rotates around (0,0,0).
#' @param throw.out.cells (Numeric) Throw out the youngest N cells in pseudotime
#' @param pseudotime (Character) Pseudotime (i.e. column name of \code{@@pseudotime}) to use for determining cells to throw out.
#'
#' @return An URD object with the coordinates of some cells in \code{@@tree$walks.force.layout} modified.
#'
#' @export
treeForceRotateCoords <- function(object, cells=NULL, seg=NULL, angle, axis=c("x", "y", "z"), around.cell=NULL, throw.out.cells=0, pseudotime=NULL) {
# If segment specified, grab the cells from that segment
if (!is.null(seg) & is.null(cells)) cells <- cellsInCluster(object, "segment", segChildrenAll(object, seg, include.self=T))
# If cells not specified, do all cells in the layout.
if (is.null(cells)) cells <- rownames(object@tree$walks.force.layout)
# Make sure all cells specified are actually in the force-directed layout
cells <- intersect(cells, rownames(object@tree$walks.force.layout))
# If throw.out.cells is on, by removing n cells with the youngest pseudotime.
if ((throw.out.cells > 0)) {
if (is.null(pseudotime)) stop ("To use throw.out.cells, must specify pseudotime.")
cells <- setdiff(cells, cells[order(object@pseudotime[cells,pseudotime])[1:throw.out.cells]])
}
# Grab the appropriate chunk of the layout.
m <- as.matrix(object@tree$walks.force.layout[cells, c("x","y","telescope.pt")])
# If a cell is specified to rotate around, translate coords so that cell is at 0, 0, 0
if (!is.null(around.cell)) {
if (is.numeric(around.cell)) {
#around.cell <- rownames(m)[order(m[,3])[around.cell]]
around.cell <- cells[order(object@pseudotime[cells,pseudotime])[around.cell]]
}
around.coords <- m[around.cell,]
m <- translateCoords3d(m=m, x=-1*around.coords[1], y=-1*around.coords[2], z=-1*around.coords[3])
}
# Do the rotation
mprime <- rotateCoords3d(m=m, a=angle, axis=axis)
# If you were rotating around an arbitrary cell, put them back
if (!is.null(around.cell)) {
mprime <- translateCoords3d(m=mprime, x=around.coords[1], y=around.coords[2], z=around.coords[3])
}
object@tree$walks.force.layout[cells, c("x", "y", "telescope.pt")] <- mprime
return(object)
}
#' Translate points in force-directed layout
#'
#' This function can be used in addition to \code{\link{treeForceRotateCoords}}
#' in order to fine-tune the presentation of a force-directed layout. For instance,
#' this can be used to achieve an improved 2D visualization by modifying overlapping
#' portions of the layout.
#'
#' @param object An URD object
#' @param cells (Character vector) Cells to modify (Default \code{NULL} is all cells in the force-directed layout)
#' @param seg (Character) Instead of specifying cells, just grab all cells from this segment and downstream. Ignored if \code{cells} is specified.
#' @param x (Numeric) Distance to move cells along x-axis
#' @param y (Numeric) Distance to move cells along y-axis
#' @param z (Numeric) Distance to move cells along z-axis
#'
#' @return An URD object with the coordinates of some cells in \code{@@tree$walks.force.layout} modified.
#'
#' @export
treeForceTranslateCoords <- function(object, cells=NULL, seg=NULL, x, y, z) {
if (!is.null(seg) & is.null(cells)) cells <- cellsInCluster(object, "segment", segChildrenAll(object, seg, include.self=T))
if (is.null(cells)) cells <- rownames(object@tree$walks.force.layout)
cells <- intersect(cells, rownames(object@tree$walks.force.layout))
m <- as.matrix(object@tree$walks.force.layout[cells, c("x","y","telescope.pt")])
mprime <- translateCoords3d(m=m, x=x, y=y, z=z)
object@tree$walks.force.layout[cells, c("x", "y", "telescope.pt")] <- mprime
return(object)
}
#' Stretch points in force-directed layout
#'
#' This function can be used in addition to \code{\link{treeForceRotateCoords}}
#' and \code{\link{treeForceTranslateCoords}}
#' in order to fine-tune the presentation of a force-directed layout. For instance,
#' this can be used to achieve an improved 2D visualization by shortening or lengthening
#' a branch of the tree if branches have dramatically different lengths in pseudotime.
#' (This occurs, for instance, in data sets generated from a single timepoint in
#' a homeostatic animal instead of an explicit timecourse.) This transformation
#' will only occur in the force-directed layout, not the dendrogram tree layout.
#'
#' @param object An URD object
#' @param factor (Numeric) Factor to lengthen (>1) or shrink (<1) the branch by.
#' @param cells (Character vector) Cells to modify (Default \code{NULL} is all cells in the force-directed layout)
#' @param seg (Character) Instead of specifying cells, just grab all cells from this segment and downstream. Ignored if \code{cells} is specified.
#' @param throw.out.cells (Numeric) Throw out the youngest N cells in pseudotime (i.e. don't stretch those)
#'
#' @return An URD object with the coordinates of some cells in \code{@@tree$walks.force.layout} modified.
#'
#' @export
treeForceStretchCoords <- function(object, factor, cells=NULL, seg=NULL, throw.out.cells=0) {
if (!is.null(seg) & is.null(cells)) cells <- cellsInCluster(object, "segment", segChildrenAll(object, seg, include.self=T))
if (is.null(cells)) cells <- rownames(object@tree$walks.force.layout)
cells <- intersect(cells, rownames(object@tree$walks.force.layout))
if ((throw.out.cells > 0)) {
#if (is.null(pseudotime)) stop ("To use throw.out.cells, must specify pseudotime.")
cells <- setdiff(cells, cells[order(object@tree$pseudotime[cells])[1:throw.out.cells]])
#cells <- setdiff(cells, cells[order(object@pseudotime[cells,pseudotime])[1:throw.out.cells]])
}
m <- object@tree$walks.force.layout[cells, c("x","y","telescope.pt")]
m.min <- min(m$telescope.pt)
m$telescope.pt <- ((m$telescope.pt - m.min) * factor) + m.min
object@tree$walks.force.layout[cells, c("x", "y", "telescope.pt")] <- m
return(object)
}
treeForceStretchCoords <- function(object, factor, cells=NULL, seg=NULL, stretch.from=20, throw.out.cells=0, ramp=0.1) {
if (!is.null(seg) & is.null(cells)) cells <- cellsInCluster(object, "segment", segChildrenAll(object, seg, include.self=T))
if (is.null(cells)) cells <- rownames(object@tree$walks.force.layout)
cells <- intersect(cells, rownames(object@tree$walks.force.layout))
if ((throw.out.cells > 0)) {
#if (is.null(pseudotime)) stop ("To use throw.out.cells, must specify pseudotime.")
cells <- setdiff(cells, cells[order(object@tree$pseudotime[cells])[1:throw.out.cells]])
#cells <- setdiff(cells, cells[order(object@pseudotime[cells,pseudotime])[1:throw.out.cells]])
}
m <- object@tree$walks.force.layout[cells, c("x","y","telescope.pt")]
m.min <- sort(m$telescope.pt)[stretch.from]
m.max <- max(m$telescope.pt)
m.fac <- (m$telescope.pt - m.min) / (m.max - m.min)
m.fac <- pmax(m.fac, rep(0, length(m.fac)))
m.fac <- pmin(m.fac / ramp, rep(1, length(m.fac))) * factor
m$telescope.pt <- ((m$telescope.pt - m.min) * m.fac) + m.min
object@tree$walks.force.layout[cells, c("x", "y", "telescope.pt")] <- m
return(object)
}
farrellja/URD documentation built on June 17, 2020, 4:48 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions treeForceStretchCoords treeForceStretchCoords treeForceTranslateCoords treeForceRotateCoords translateCoords3d rotateCoords3d plotTreeForceStore3DView fdlDensity treeForcePositionLabels plotTreeForceDual plotTreeForce plotTreeForce2D treeForceDirectedLayout
Documented in fdlDensity plotTreeForce plotTreeForce2D plotTreeForceDual plotTreeForceStore3DView rotateCoords3d translateCoords3d treeForceDirectedLayout treeForcePositionLabels treeForceRotateCoords treeForceStretchCoords treeForceTranslateCoords
#' Generate force-directed layout using tip-walk data.
#'
#' This function generates a k-nearest neighbor network for cells, based on their visitation
#' frequency by the biased random walks from different tips (and optionally pseudotime), and
#' uses it as input into a force-directed layout (powered by igraph). The force-directed
#' layout is generated in 2 dimensions, and pseudotime is used as a third dimension.
#'
#' Several settings adjust the k-nearest neighbor network to assist in the layout. Outlier
#' cells (that are only closely connected to a single other cells) or outlier edges (with
#' unusually long distances) can be eliminated, though these parameters are disabled by
#' default. However, the dendrogram structure recovered by URD is used by default to
#' refine the k-nearest neighbor network, breaking links between cells that are distant
#' in the dendrogram. This emphasizes producing a parseable layout at the expense of
#' representing rare transitions in the data.
#'
#' Parameters that greatly affect the quality of the layout are the number of nearest
#' neighbors used in the graph (\code{num.nn}) and the aggresiveness of the
#' \code{cut.unconnected.segments} parameter. Additionally, the layout is largely
#' reproducible if given the same starting conditions, but minor changes can affect
#' it quite a bit, so it is sometimes worth varying \code{num.nn} by 1 across a small range
#' (i.e. try 5 layouts with \code{num.nn=120-125}) and choosing the best one. Finally,
#' sometimes different branches can overlap after they have separated, and for a final
#' publication-ready layout, it can make sense to hand-tune the layout some. This can
#' be done using the \code{\link{treeForceRotateCoords}} and \code{\link{treeForceTranslateCoords}}
#' functions to move or rotate sections of the tree. (For instance, in the zebrafish
#' layout, we used these functions to move the completely disconnected EVL and PGC cells
#' into place and to increase the angular distance at the first branchpoint in the blastoderm
#' so that the more fine-grained branches from the different germ layers didn't overlap.)
#'
#' Thanks for Dorde Relic for debugging and creating the \code{remove.duplicate.cells} parameter and Simon Cai for debugging that led to the \code{cell.minimum.walks} parameter.
#'
#' @importFrom RANN nn2
#' @importFrom stats quantile
#' @importFrom reshape2 melt
#' @importFrom igraph graph_from_data_frame layout_with_fr layout_with_drl layout_with_kk vcount V distances is_connected decompose
#'
#' @param object An URD object
#' @param num.nn (Numeric) Number of nearest neighbors to use. (\code{NULL} will use the square root of the number of cells as default.)
#' @param method (Character: "fr", "drl", or "kk") Which force-directed layout algorithm to use (\link[igraph:layout_with_fr]{Fruchterman-Reingold}, \link[igraph:layout_with_drl]{DrL}, or \link[igraph:layout_with_kk]{Kamada-Kawai})
#' @param cells.to.do (Character vector) Cells to use in the layout (default \code{NULL} is all cells in the tree.)
#' @param cell.minimum.walks (Numeric) Minimum number of times a cell must have been visited by random walks in order to be included in the force-directed layout (Cells that have been visited only a few times are more likely to be assigned incorrectly or poorly.)
#' @param cut.outlier.cells (Numeric) If desired, omit cells with unusual second nearest neighbor distances (i.e. those that are likely outliers only well connected to one other cell). Parameter is given as a factor of the interquartile range calculated across all cells' distance to their second nearest neighbor. (Default is not to omit any cells.)
#' @param cut.outlier.edges (Numeric) If desired, cut edges in the nearest neighbor graph with unusually long distances. Parameters is given as a factor of the interquartile range calculated across all edges in the graph. (Default is not to cut any edges based on their length.)
#' @param max.pseudotime.diff (Numeric) If desired, cut edges in the nearest neighbor graph between cells with longer difference in pseudotime. (Default is not to cut any edges based on their pseudotime.)
#' @param cut.unconnected.segments (Numeric) Cut connections in the nearest-neighbor graph to cells that are
#' more segments away in the dendrogram structure. For instance, the most aggressive setting (1) will only
#' maintain that are at most 1 segment away (so it will maintain connections within a segment and to that
#' segment's parent or children. Higher values permit more distant connections. The default value is 2,
#' which permits connections up to two segments away (i.e. connections within a segment, to that segment's
#' parent, grandparent, children, grandchildren, and siblings.) \code{NA} or \code{NULL} disable this setting
#' and permit all connections.
#' @param min.final.neighbors (Numeric) After trimming outlier and unconnected connections in the nearest
#' neighbor graph, remove any cells that remain connected to fewer than this many other cells.
#' @param remove.duplicate.cells (Logical) If cells have exactly the same visitation during the random walks from each tip, this can create problems in the force-directed layout where cells are pushed to the outside of the layout. To avoid this if \code{remove.duplicate.cells=T}, only one cell from each group of cells with duplicated coordinates will be used in the layout.
#' @param tips (Character vector) Tips for which walk visitation data should be used in the construction of the nearest neighbor graph. (Default is all tips )
#' @param coords (Matrix: Cells as rows, 2 columns) Starting coordinates for the force directed layout. Default (\code{"auto"}) takes them from the cell layout of the dendrogram.
#' @param start.temp (Numeric) Starting temperature for the force-directed layout (if \code{method="fr"}), which controls how much cells can move in the initial iterations of the algorithm. Default (\code{NULL}) is the square root of the number of cells.
#' @param n.iter (Numeric or \code{NULL}) Sets the number (or maximum number) of iterations used in the force-directed layout; if \code{NULL}, observes defaults suggested by igraph.
#' @param density.neighbors (Numeric) Distance to this nearest neighbor (default is 10th nearest neighbor) is used as a proxy for local density in the force-directed layout. This is used by \code{\link{plotTreeForce}} if \code{density.alpha=T} for increasing transparency in more high density regions of the layout. This can be adjusted after generating the layout by re-running the \code{\link{fdlDensity}} function.
#' @param plot.outlier.cuts (Logical) If \code{cut.outlier.cells=T} or \code{cut.outlier.edges=T}, this displays the
#' @param verbose (Logical) Print progress and time stamps?
#'
#' @examples
#' object.built <- treeForceDirectedLayout(object.built, num.nn=120, pseudotime="pseudotime", method = "fr", dim = 2, cells.to.do = robustly.visited.cells, tips=final.tips, cut.unconnected.segments = 2, min.final.neighbors=4, verbose=T)
#'
#' @export
treeForceDirectedLayout <- function(object, num.nn=NULL, method=c("fr", "drl", "kk"), cells.to.do=NULL, cell.minimum.walks=1, cut.outlier.cells=NULL, cut.outlier.edges=NULL, max.pseudotime.diff=NULL, cut.unconnected.segments=2, min.final.neighbors=2, remove.duplicate.cells=T, tips=object@tree$tips, coords="auto", start.temp=NULL, n.iter=NULL, density.neighbors=10, plot.outlier.cuts=F, verbose=F) {
# Params
if (length(method) > 1) method <- method[1]
if (is.null(cells.to.do)) cells.to.do <- rownames(object@diff.data)
starting.cells <- length(cells.to.do)
if (is.null(num.nn)) num.nn <- ceiling(sqrt(length(cells.to.do)))
if (is.null(tips)) tips <- object@tree$tips
if (class(coords)=="character" && coords=="auto") {
coords <- as.matrix(object@tree$cell.layout[cells.to.do,c("x","y")])
} else if (!is.null(coords) && class(coords) != "matrix") {
stop("coords must either be 'auto', NULL, or a matrix.")
}
dim=2
if (verbose) print(paste(Sys.time(), ": Starting with parameters", method, num.nn, "NN", dim, "D", length(cells.to.do), "cells"))
# Get rid of any cells that don't have a pseudotime (or, if cut.unconnected.segments is turned on, any cells that aren't part of the tree)
cells.no.pseudotime <- cells.to.do[which(is.na(object@tree$pseudotime[cells.to.do]))]
if (!is.na(cut.unconnected.segments) & !is.null(cut.unconnected.segments) & cut.unconnected.segments > 0) {
cells.not.in.tree <- setdiff(cells.to.do, unlist(object@tree$cells.in.segment))
if (verbose) print(paste0("Removing ", length(unique(c(cells.no.pseudotime, cells.not.in.tree))), " cells that are not assigned a pseudotime or a segment in the tree."))
cells.to.do <- setdiff(cells.to.do, c(cells.no.pseudotime, cells.not.in.tree))
} else {
if (verbose) print(paste0("Removing ", length(cells.no.pseudotime), " cells that are not assigned a pseudotime."))
cells.to.do <- setdiff(cells.to.do, cells.no.pseudotime)
}
# Get and normalize walk data, then add pseudotime (unless NULL in which case don't use)
if (verbose) print(paste0(Sys.time(), ": Preparing walk data."))
walk.data <- object@diff.data[cells.to.do,paste0("visitfreq.raw.", object@tree$tips)]
walk.data$pseudotime <- object@tree$pseudotime[cells.to.do]
walk.total <- apply(walk.data, 1, sum)
# If some cells weren't walked enough
if (any(walk.total < cell.minimum.walks)) {
if (verbose) print(paste0("Removing ", length(which(walk.total < cell.minimum.walks)), " cells that were visited fewer than ", cell.minimum.walks, " times by random walks."))
# Remove those cells
walk.data <- walk.data[which(walk.total >= cell.minimum.walks),]
# And recalculate walk totals
walk.total <- apply(walk.data, 1, sum)
}
walk.data <- sweep(walk.data, 1, walk.total, "/")
walk.data <- as.matrix(walk.data)
duped <- which(duplicated(walk.data))
if (length(duped) > 0 && remove.duplicate.cells) {
warning(length(duped), " cells have duplicate random walk coordinates and are being removed from the layout.")
walk.data <- unique(walk.data)
} else if (length(duped) > 0) {
warning(length(duped), " cells have duplicate random walk coordinates. This may cause a problem in the layout. If so, set remove.duplicate.cells=T")
}
# Calculate nearest neighbor graph
if (verbose) print(paste0(Sys.time(), ": Calculating nearest neighbor graph."))
walk.nn <- RANN::nn2(data=walk.data,query=walk.data,k=max(num.nn)+1, treetype = "kd", searchtype="priority")
rownames(walk.nn$nn.idx) <- rownames(walk.data)
walk.nn$nn.idx <- walk.nn$nn.idx[,2:(num.nn+1)]
rownames(walk.nn$nn.dists) <- rownames(walk.data)
walk.nn$nn.dists <- walk.nn$nn.dists[,2:(num.nn+1)]
walk.nn$nn.label <- apply(walk.nn$nn.idx, 2, function(y) rownames(walk.data)[y])
rownames(walk.nn$nn.label) <- rownames(walk.data)
# Cull outlier cells with unusual second nearest neighbor distances
if (!is.null(cut.outlier.cells)) {
q <- quantile(walk.nn$nn.dists[,2])
outer.fence <- q[4] + cut.outlier.cells*(q[4]-q[1])
if (plot.outlier.cuts) {
hist(walk.nn$nn.dists[,2], breaks=100, main="Outlier Removal: Second Nearest Neighbor", xlab="Distance to 2nd Nearest Neighbor")
abline(v=outer.fence, col='red')
}
if (verbose) print(paste0(Sys.time(), ": Removing ", round(length(which(walk.nn$nn.dists[,2] > outer.fence)) / starting.cells, digits=2), "% of cells as outliers."))
cells.keep <- which(walk.nn$nn.dists[,2] <= outer.fence)
} else {
cells.keep <- 1:dim(walk.data)[1]
}
# Turn nearest neighbors into edge list
if (verbose) print(paste0(Sys.time(), ": Preparing edge list."))
edges <- reshape2::melt(walk.nn$nn.label[cells.keep,], stringsAsFactors=F)
dists <- reshape2::melt(walk.nn$nn.dists[cells.keep,])
edges <- edges[,c(1,3)]
names(edges) <- c("V1", "V2")
edges$dists <- dists$value
edges$V1 <- as.character(edges$V1)
edges$V2 <- as.character(edges$V2)
starting.edges <- dim(edges)[1]
# Cut connections with unusually long distances
if (!is.null(cut.outlier.edges)) {
q <- quantile(edges$dists)
outer.fence <- q[4] + cut.outlier.edges*(q[4]-q[1])
if (verbose) print(paste0(Sys.time(), ": Removing ", round(length(which(edges$dists > outer.fence)) / starting.edges * 100, digits=3), "% of edges as outliers."))
if (plot.outlier.cuts) {
hist(edges$dists, breaks=100, main="Outlier Removal: Edge distances", xlab="Edge distance")
abline(v=outer.fence, col='red')
}
edges <- edges[edges$dists <= outer.fence,]
}
# Cut connections with overly large pseudotime differences
if (!is.null(max.pseudotime.diff)) {
edges$pt1 <- object@tree$pseudotime[edges$V1]
edges$pt2 <- object@tree$pseudotime[edges$V2]
edges$dpt <- abs(edges$pt1 - edges$pt2)
if (verbose) print(paste0(Sys.time(), ": Removing ", round(length(which(edges$dpt > max.pseudotime.diff)) / starting.edges * 100, digits=3), "% of edges with too large pseudotime difference."))
edges <- edges[which(edges$dpt <= max.pseudotime.diff),]
}
# Cut edges between overly distant (e.g. non-connected) segments
if (!is.na(cut.unconnected.segments) & !is.null(cut.unconnected.segments) & cut.unconnected.segments > 0) {
# Determine which segments each connected cell belongs to
edges$seg1 <- object@diff.data[edges$V1, "segment"]
edges$seg2 <- object@diff.data[edges$V2, "segment"]
# Determine distance in the tree between each segment
seg.dist <- igraph::distances(object@tree$tree.igraph, mode="all")
edges$seg.dist <- apply(edges[,c("seg1","seg2")], 1, function(segs) seg.dist[segs[1],segs[2]])
if (verbose) print(paste0(Sys.time(), ": Removing ", round(length(which(edges$seg.dist > cut.unconnected.segments)) / starting.edges * 100, digits=2), "% of edges that are between segments with distance > ", cut.unconnected.segments))
edges <- edges[edges$seg.dist <= cut.unconnected.segments,]
}
# Trim cells that are no longer well connected to the graph
if (verbose) print(paste0(Sys.time(), ": Trimming cells that are no longer well connected."))
connections.remaining <- table(edges[edges$dists > 0,"V1"])
cells.with.enough.connections <- names(which(connections.remaining >= min.final.neighbors))
cells.without.enough.connections <- length(unique(edges$V1)) - length(cells.with.enough.connections)
while(cells.without.enough.connections > 0) {
edges <- edges[which(edges$V1 %in% cells.with.enough.connections & edges$V2 %in% cells.with.enough.connections),]
connections.remaining <- table(edges[edges$dists > 0,"V1"])
cells.with.enough.connections <- names(which(connections.remaining >= min.final.neighbors))
cells.without.enough.connections <- length(unique(edges$V1)) - length(cells.with.enough.connections)
}
if (verbose) print(paste0(Sys.time(), ": ", round(length(cells.with.enough.connections)/starting.cells*100, digits=2), "% of starting cells preserved."))
# Create igraph representation
if (verbose) print(paste0(Sys.time(), ": Preparing igraph object."))
if (method=="kk") {
# kk requires that weights be higher for points that should be far apart -- use distance directly
edges$weight <- edges$dists
} else {
# everything else places points with high weight nearby -- use inverse of distance as graph weight
edges$weight <- 1/edges$dists
}
object@tree$walks.force.edges <- edges
igraph.walk.weights <- igraph::graph_from_data_frame(edges, directed=F)
if (!is.null(coords)) {
cells.remain <- unique(c(edges$V1, edges$V2))
coords.remain <- rownames(coords) %in% cells.remain
coords <- coords[coords.remain,]
}
# DrL force-directed layout only works on connected graphs, so check for graph connectivity
if (method == "drl" && !igraph::is_connected(igraph.walk.weights)) {
igraph.walk.weights.decomposed <- decompose(igraph.walk.weights)
graph.keep <- which.max(unlist(lapply(igraph.walk.weights.decomposed, function(x) length(V(x)))))
warning("DrL method can only operate on fully connected graphs. Using more nearest neighbors will increase graph connectivity. For now, using largest subgraph: ", length(V(igraph.walk.weights.decomposed[[graph.keep]])), " of ", length(V(igraph.walk.weights)), " cells.")
igraph.walk.weights <- igraph.walk.weights.decomposed[[graph.keep]]
coords <- coords[rownames(coords) %in% names(V(igraph.walk.weights)),]
cells.with.enough.connections <- intersect(cells.with.enough.connections, names(V(igraph.walk.weights)))
}
# Define start temp, if set to default
if (is.null(start.temp)) start.temp <- sqrt(igraph::vcount(igraph.walk.weights))
# Define number of iterations to do
if (is.null(n.iter)) {
if (method == "fr") n.iter = 500
if (method == "kk") n.iter = 50*igraph::vcount(igraph.walk.weights)
}
# Do force-directed layout
if (verbose) print(paste0(Sys.time(), ": Doing force-directed layout."))
if (method=="fr") igraph.walk.layout <- igraph::layout_with_fr(igraph.walk.weights, dim=dim, coords=coords, start.temp=start.temp, niter=n.iter)
if (method=="drl") igraph.walk.layout <- igraph::layout_with_drl(igraph.walk.weights, dim=dim, options=list(edge.cut=0))
if (method=="kk") igraph.walk.layout <- igraph::layout_with_kk(igraph.walk.weights, dim=dim, coords=coords, maxiter=n.iter)
# Store the layout
if (dim==2) {
# Calculate weighted pseudotime of cells' neighbors for tree layout.
if (verbose) print(paste0(Sys.time(), ": Calculating Z."))
neighbor.pt <- unlist(lapply(cells.with.enough.connections, function(i) {
# Deal with points with 0 distance, by setting their weight to the sum of all other weights.
weights <- 1/walk.nn$nn.dists[i,]
weights[!is.finite(weights)] <- sum(weights[is.finite(weights)])
weighted.mean(x=object@tree$pseudotime[walk.nn$nn.label[i,]], w=weights)
}))
names(neighbor.pt) <- cells.with.enough.connections
# Store the layout
object@tree$walks.force.layout <- as.data.frame(igraph.walk.layout, stringsAsFactors=F)
names(object@tree$walks.force.layout) <- c("x","y")
rownames(object@tree$walks.force.layout) <- igraph::V(igraph.walk.weights)$name
# Remove points with NA values for force-directed layout
fdl.failed <- length(which(!(complete.cases(object@tree$walks.force.layout))))
if (fdl.failed > 0) {
warning(fdl.failed, " cells were not assigned coordinates during the force-directed layout.")
object@tree$walks.force.layout <- object@tree$walks.force.layout[complete.cases(object@tree$walks.force.layout),]
}
# Normalize neighbor.pt range
neighbor.pt.factor <- mean(apply(object@tree$walks.force.layout[,c("x","y")], 2, max)) / max(neighbor.pt, na.rm=T)
object@tree$walks.force.layout$telescope.pt <- neighbor.pt[rownames(object@tree$walks.force.layout)] * neighbor.pt.factor
# Calculate local density for adjusting alpha
if (verbose) print(paste0(Sys.time(), ": Calculating local density."))
object <- fdlDensity(object, neighbor=density.neighbors)
} else if (dim==3) {
object@tree$walks.force.layout.3d <- as.data.frame(igraph.walk.layout, stringsAsFactors=F)
names(object@tree$walks.force.layout.3d) <- c("x","y","z")
rownames(object@tree$walks.force.layout.3d) <- igraph::V(igraph.walk.weights)$name
}
# Position labels if segments have been named.
if (!is.null(object@tree$segment.names)) {
object@tree$walks.force.labels <- treeForcePositionLabels(object)
}
# Return the object
if (verbose) print(paste0(Sys.time(), ": Finished."))
return(object)
}
#' Plot tree force-directed layout in 2D
#'
#' @import ggplot2
#'
#' @param object An URD object
#'
#' @export
plotTreeForce2D <- function(object, label=NULL, label.type="search", title=label, show.points=T, point.alpha=1, point.size=1, show.neighbors=F, neighbors.max=10000, colors=NULL) {
# Load colors if not specified
if (is.null(colors)) colors <- defaultURDContinuousColors()
gg.data <- object@tree$walks.force.layout
if (!is.null(label)) color.data <- data.for.plot(object = object, label = label, label.type = label.type, as.color = F, as.discrete.list=T, cells.use = rownames(gg.data))
gg.data$expression <- color.data$data
gg <- ggplot() + theme_bw() + labs(x="", y="", color="", title=title)
if (show.neighbors) {
edge.data <- object@tree$walks.force.edges
if (dim(edge.data)[1] > neighbors.max) edge.data <- edge.data[sample(1:dim(edge.data)[1], neighbors.max, replace=F),]
edge.data[,c("x1","y1")] <- gg.data[edge.data$V1,c("x","y")]
edge.data[,c("x2","y2")] <- gg.data[edge.data$V2,c("x","y")]
if (show.points) gg <- gg + geom_segment(data=edge.data, aes(x=x1, y=y1, xend=x2, yend=y2), color='black', alpha=0.02)
if (!show.points) gg <- gg + geom_segment(data=edge.data, aes(x=x1, y=y1, xend=x2, yend=y2, color=weight), alpha=0.04)
}
if (show.points) gg <- gg + geom_point(data=gg.data, aes(x=x, y=y, color=expression), size=point.size, alpha=point.alpha)
if (!color.data$discrete) gg <- gg + scale_color_gradientn(colors = colors)
return(gg)
}
#' Plot force-directed layout of tree
#'
#' This plots data on the force-directed layout of an URD tree, generated by
#' \code{\link{treeForceDirectedLayout}}. Once the tree is rotated into a
#' perspective that is appealing, \code{\link{plotTreeForceStore3DView}} can be
#' used to store that view in the URD object for additional plots.
#'
#' The transparency of points is typically adjusted in two ways: (1) Points with
#' low expression are made more transparent, which is controlled through the
#' \code{fade.below} and \code{alpha.fade} parameters. (2) Points in higher
#' density regions of the plot are made more transparent, which is controlled by
#' the \code{density.alpha} parameter (and the local density is determined by
#' the \code{density.neighbors} parameter of \code{\link{treeForceDirectedLayout}},
#' and can be adjusted after the fact by re-running \code{\link{fdlDensity}}).
#'
#' Plots can be easily saved using \code{\link{rgl:rgl.snapshot}} and plots
#' can be closed with \code{\link{rgl:rg.close}} if plotting many genes
#' using a loop.
#'
#' @param object An URD object
#' @param label (Character) Data to use for coloring tree (see \link{data.for.plot})
#' @param label.type (Character) Where to find data for label (default \code{"search"} auto-detects, see \link{data.for.plot})
#' @param view (Character) Name of view to use (i.e. the name of an entry in \code{@@tree$force.view.list}) created by \link{plotTreeForceStore3DView}. Default is last view created.
#' @param alpha (Numeric) Maximum transparency of points
#' @param alpha.fade (Numeric) Minimum transparency of points
#' @param size (Numeric) Size of points in the plot
#' @param title (Character) Title to add to the plot. (This is sensitive to resizing the window after plotting, but if a view is stored, the window will be resized before the title is added, so it will be acceptable resolution for figures.)
#' @param title.cex (Numeric) Adjust the title font size
#' @param title.line (Numeric) Adjust the position of the title. Positive numbers move the title upward.
#' @param label.tips (Logical) Should text identifying the tips of the tree be placed in the force directed layout? Defaults to TRUE if \code{@@tree$segment.names} has been set.
#' @param use.short.names (Logical) Should short names from \code{@@tree$segment.names.short} be used? Defaults to TRUE if those values have been set.
#' @param seg.show (Character vector) Segments of the tree to put on the plot (Default \code{NULL} is all segments)
#' @param cells.show (Character vector) Cells of the tree to show (Default \code{NULL} is all cells)
#' @param fade.below (Numeric) If desired, transparency can be lowered for points with low expression values (thereby highlighting the portions of the tree with expression). This value is the portion of the range of expression values to fade. (The default, 2/9, fades the bottom two-ninths of the expression values.) A fully 'faded' point will have alpha of \code{alpha.fade}.
#' @param density.alpha (Logical) Should points with short distance to nearest neighbors (i.e. in denser regions of the plot) be more transparent?
#' @param label.spacing (Numeric) How far should labels be from the detected tips
#' @param text.cex (Numeric) Size of the label text
#' @param colors (Character vector) Vector of colors to use if plotting continuous data
#' @param discrete.colors (Character vector) Vector of colors to use if plotting discrete data
#' @param color.limits (Numeric vector, length 2) Minimum and maximum values for color scale. Default \code{NULL} auto-detects.
#' @param symmetric.color.scale (Logical) Should the color scale be symmetric and centered around 0? (Default \code{NULL} is \code{FALSE} if all values are positive, and \code{TRUE} if both positive and negative values are present.)
#'
#' @return Nothing. Produces a plot using the \code{rgl} package, displayed in an X11 window.
#'
#' @examples
#' # Define the colors used below
#' stage.colors <- c("#CCCCCC", RColorBrewer::brewer.pal(9, "Set1")[9], RColorBrewer::brewer.pal(12, "Paired")[c(9,10,7,8,5,6,3,4,1,2)])
#' fire.with.grey <- c("#CECECE", "#DDC998", RColorBrewer::brewer.pal(9, "YlOrRd")[3:9])
#' branch.colors <- c("#CECECE", "#E6298B")
#'
#' # Plot a discrete characteristic like stage or clustering results
#' plotTreeForce(object, "stage.nice", title="STAGE", title.line=1, discrete.colors = stage.colors, alpha=0.4)
#'
#' # Plot gene expression
#' plotTreeForce(object, "NOTO", title="NOTO expression", title.line=1, colors = fire.with.grey)
#'
#' Define a 'clustering' based on whether or not cells are in a given trajectory then plot it
#' object <- groupFromCells(object, group.id="lineage_Tailbud", cells=cellsAlongLineage(object, "Tailbud", remove.root=F))
#' plotTreeForce(object, "lineage_Tailbud", title="Tailbud", title.line=1, discrete.colors=branch.colors, alpha=0.4)
#'
#' @export
plotTreeForce <- function(object, label, label.type="search", view="default", alpha=0.8, alpha.fade=0.1, size=5, title=NULL, title.cex=3, title.line=0, label.tips=(!is.null(object@tree$segment.names) | !is.null(object@tree$segment.names.short)), use.short.names=!is.null(object@tree$segment.names.short), seg.show=NULL, cells.show=NULL, fade.below=(2/9), density.alpha=T, label.spacing=5, text.cex=0.8, colors=NULL, discrete.colors=NULL, color.limits=NULL, symmetric.color.scale=NULL) {
if (requireNamespace("rgl", quietly = TRUE)) {
# Figure out which view to use
if (!is.null(view)) {
if (view=="default") {
if ("force.view.default" %in% names(object@tree)) view <- object@tree$force.view.default
else view <- NULL
}
view <- object@tree$force.view.list[[view]]
}
# Get layout data
gg.data <- object@tree$walks.force.layout
gg.data$segment <- object@diff.data[rownames(gg.data),"segment"]
# Get expression data, use to determine parameters of color scale, and get default colors
expression.data <- data.for.plot(object=object, label=label, label.type=label.type, as.color=F, cells.use=rownames(gg.data), as.discrete.list=T)
if (!expression.data$discrete && is.null(symmetric.color.scale)) symmetric.color.scale <- min(expression.data$data) < 0
if (!expression.data$discrete && is.null(color.limits)) {
if (symmetric.color.scale) {
color.sv <- max(abs(expression.data$data))
color.limits <- c(-1*color.sv, color.sv)
} else {
color.limits <- c(0, max(expression.data$data))
}
}
if (is.null(colors) && !expression.data$discrete) colors <- defaultURDContinuousColors(symmetric = symmetric.color.scale)
# Set color, alpha, and point size
if (!is.null(label)) {
color.data <- data.for.plot(object = object, label = label, label.type = label.type, as.color = T, as.discrete.list=T, cells.use = rownames(gg.data), continuous.colors=colors, colors.use = discrete.colors, continuous.color.limits = color.limits)
gg.data$expression <- color.data$data
}
gg.data$alpha <- alpha
gg.data$size <- size
# Focus on expression
# For this, grey out any thing below focus, set alpha low for those, and increase progressively until 2*focus.
if (fade.below > 0 & !color.data$discrete) {
# Figure out ranges
expression.range <- range(abs(expression.data$data))
fade.range <- diff(expression.range) * fade.below + expression.range[1]
# Adjust gg.data
fade <- (fade.range-abs(expression.data$data)+expression.range[1])/fade.range
fade[fade < 0] <- 0
gg.data$alpha <- alpha - ((alpha-alpha.fade) * fade)
}
# Open 3D view
if (!is.null(view)) {
rgl::open3d(
zoom=view$rgl.setting$zoom,
scale=view$rgl.setting$scale,
userMatrix=view$rgl.setting$userMatrix,
windowRect=view$rgl.setting$windowRect
)
} else {
rgl::open3d()
}
# Adjust transparency for local density
if (density.alpha) {
dr <- range(object@tree$walks.force.layout$n.dist)
n.dist <- object@tree$walks.force.layout$n.dist / dr[2]
density.adj <- sqrt(n.dist)
density.adj <- density.adj / median(density.adj)
gg.data$alpha <- gg.data$alpha * density.adj
}
# Convert alpha to 0 on cells not to show (so the plot size doesn't get changed)
if (!is.null(seg.show) & is.null(cells.show)) {
segs.show <- segChildrenAll(object, seg.show, include.self=T)
cells.show <- rownames(gg.data)[which(gg.data$segment %in% segs.show)]
}
if (!is.null(cells.show)) {
cells.hide <- setdiff(rownames(gg.data), cells.show)
gg.data[cells.hide,"alpha"] <- 0
}
# Add the points to the plot
rgl::points3d(x=gg.data$x, y=gg.data$y, z=gg.data$telescope.pt, col=gg.data$expression, alpha=gg.data$alpha, size=size)
# Add title to plot
if (!is.null(title)) {
Sys.sleep(0.1)
rgl::bgplot3d({plot.new(); title(main=title, line=title.line, cex.main=title.cex)})
Sys.sleep(0.1)
}
# Add labels to tips
if (label.tips) {
# Check for existing tree labels
if (!is.null(object@tree$walks.force.labels)) {
tip.labels <- object@tree$walks.force.labels
} else {
# If there aren't any, generate some.
tip.labels <- treeForcePositionLabels(object)
}
# Choose short or long labels
if (use.short.names) tip.labels$label <- tip.labels$name.short else tip.labels$label <- tip.labels$name
# Limit to labels of cells that appear in the plot
segs.to.label <- intersect(unique(gg.data[which(gg.data$alpha > 0),"segment"]), tip.labels$seg)
tip.labels <- tip.labels[which(tip.labels$seg %in% segs.to.label),]
# Add the text
rgl::text3d(x=tip.labels$x, y=tip.labels$y, z=tip.labels$z, text=tip.labels$label, adj=0.5, cex=text.cex)
}
# Add a brief pause to ensure that rendering completes before moving on to next function.
Sys.sleep(0.1)
} else {
stop("Package rgl is required for this function. To install: install.packages('rgl')\n")
}
}
#' Plot force-directed layout of tree
#'
#' This plots data on the force-directed layout of an URD tree, generated by
#' \code{\link{treeForceDirectedLayout}}, using red and green (which become yellow
#' when they overlap to display two pieces of data simultaneously. See
#' \code{\link{plotTreeForce}} for more details.
#'
#' @param object An URD object
#' @param label.red (Character) Data to use for coloring tree (red channel) (see \link{data.for.plot})
#' @param label.green (Character) Data to use for coloring tree (green channel) (see \link{data.for.plot})
#' @param label.type.red (Character) Where to find data for red channel label (default \code{"search"} auto-detects, see \link{data.for.plot})
#' @param label.type.green (Character) Where to find data for green channel label (default \code{"search"} auto-detects, see \link{data.for.plot})
#' @param view (Character) Name of view to use (i.e. the name of an entry in \code{@@tree$force.view.list}) created by \link{plotTreeForceStore3DView}. Default is last view created.
#' @param alpha (Numeric) Maximum transparency of points
#' @param alpha.fade (Numeric) Minimum transparency of points
#' @param size (Numeric) Size of points in the plot
#' @param title (Character) Title to add to the plot. (This is sensitive to resizing the window after plotting, but if a view is stored, the window will be resized before the title is added, so it will be acceptable resolution for figures.)
#' @param title.cex (Numeric) Adjust the title font size
#' @param title.line (Numeric) Adjust the position of the title. Positive numbers move the title upward.
#' @param label.tips (Logical) Should text identifying the tips of the tree be placed in the force directed layout? Defaults to TRUE if \code{@@tree$segment.names} has been set.
#' @param use.short.names (Logical) Should short names from \code{@@tree$segment.names.short} be used? Defaults to TRUE if those values have been set.
#' @param seg.show (Character vector) Segments of the tree to put on the plot (Default \code{NULL} is all segments)
#' @param cells.show (Character vector) Cells of the tree to show (Default \code{NULL} is all cells)
#' @param fade.below (Numeric) If desired, transparency can be lowered for points with low expression values (thereby highlighting the portions of the tree with expression). This value is the portion of the range of expression values to fade. (The default, 1/3, fades the bottom third of the expression values.) A fully 'faded' point will have alpha of \code{alpha.fade}.
#' @param fade.method (Character: "max", "min", or "mean") Should fade be calculated on max expression of the two labels (emphasize any cells expressing either label), min expression of the two labels (emphasize only cells co-expressing both labels), or mean expression
#' @param density.alpha (Logical) Should points with short distance to nearest neighbors (i.e. in denser regions of the plot) be more transparent?
#' @param label.spacing (Numeric) How far should labels be from the detected tips
#' @param text.cex (Numeric) Size of the label text
#'
#' @return Nothing. Produces a plot using the \code{rgl} package, displayed in an X11 window.
#'
#' @examples
#' # Focus on all cells that express either gene
#' plotTreeForceDual(object, label.red="SOX17", label.green="SOX32", title="SOX32 (green) + SOX17 (red)")
#'
#' # Focus on only the co-expressing cells by using fade.method="min"
#' plotTreeForceDual(object, label.red="SOX17", label.green="SOX32", title="SOX32 (green) + SOX17 (red)", fade.method="min")
#'
#' @export
plotTreeForceDual <- function(object, label.red, label.green, label.red.type="search", label.green.type="search", view="default", alpha=0.8, alpha.fade=0.025, size=5, title=NULL, title.cex=3, title.line=0, label.tips=(!is.null(object@tree$segment.names) | !is.null(object@tree$segment.names.short)), use.short.names=!is.null(object@tree$segment.names.short), seg.show=NULL, cells.show=NULL, fade.below=(1/3), fade.method=c("max", "min", "mean"), density.alpha=T, label.spacing=5, text.cex=0.8) {
if (requireNamespace("rgl", quietly = TRUE)) {
# Figure out which view to use
if (!is.null(view)) {
if (view=="default") {
if ("force.view.default" %in% names(object@tree)) view <- object@tree$force.view.default
else view <- NULL
}
view <- object@tree$force.view.list[[view]]
}
# Get layout data
gg.data <- object@tree$walks.force.layout
gg.data$segment <- object@diff.data[rownames(gg.data),"segment"]
# Get expression data
red.data <- data.for.plot(object=object, label=label.red, label.type=label.red.type, as.single.color = T, as.discrete.list=T, cells.use=rownames(gg.data))
green.data <- data.for.plot(object=object, label=label.green, label.type=label.green.type, as.single.color = T, as.discrete.list=T, cells.use=rownames(gg.data))
if (red.data$discrete | green.data$discrete) stop("For dual-color plots, both labels must be continuous values.")
gg.data$expression <- rgb(red.data$data, green.data$data, 0)
gg.data$alpha <- alpha
gg.data$size <- size
# Focus on expression
# For this, grey out any thing below focus, set alpha low for those, and increase progressively until 2*focus.
if (fade.below > 0) {
# Figure out value for fading
if (length(fade.method) > 1) fade.method <- fade.method[1]
if (tolower(fade.method) == "max") {
# Use max of red or green for fading
expression.data <- pmax(red.data$data, green.data$data)
} else if (tolower(fade.method) == "min") {
# Use min of red or green for fading
expression.data <- pmin(red.data$data, green.data$data)
} else if (tolower(fade.method) == "mean") {
# Use min of red or green for fading
expression.data <- apply(data.frame(red.data$data, green.data$data), 1, mean)
} else {
stop("fade.method must be max, min, or mean.")
}
# Adjust gg.data
fade <- (fade.below - expression.data)/fade.below
fade[fade < 0] <- 0
gg.data$alpha <- alpha - ((alpha-alpha.fade) * fade)
}
# Open 3D view
if (!is.null(view)) {
rgl::open3d(
zoom=view$rgl.setting$zoom,
scale=view$rgl.setting$scale,
userMatrix=view$rgl.setting$userMatrix,
windowRect=view$rgl.setting$windowRect
)
} else {
rgl::open3d()
}
# Adjust transparency for local density
if (density.alpha) {
dr <- range(object@tree$walks.force.layout$n.dist)
n.dist <- object@tree$walks.force.layout$n.dist / dr[2]
density.adj <- sqrt(n.dist)
density.adj <- density.adj / median(density.adj)
gg.data$alpha <- gg.data$alpha * density.adj
}
# Convert alpha to 0 on cells not to show (so the plot size doesn't get changed)
if (!is.null(seg.show) & is.null(cells.show)) {
segs.show <- segChildrenAll(object, seg.show, include.self=T)
cells.show <- rownames(gg.data)[which(gg.data$segment %in% segs.show)]
}
if (!is.null(cells.show)) {
cells.hide <- setdiff(rownames(gg.data), cells.show)
gg.data[cells.hide,"alpha"] <- 0
}
# Add the points to the plot
rgl::points3d(x=gg.data$x, y=gg.data$y, z=gg.data$telescope.pt, col=gg.data$expression, alpha=gg.data$alpha, size=size)
# Add title to plot
if (!is.null(title)) {
Sys.sleep(0.1)
rgl::bgplot3d({plot.new(); title(main=title, line=title.line, cex.main=title.cex)})
Sys.sleep(0.1)
}
# Add labels to tips
if (label.tips) {
# Check for existing tree labels
if (!is.null(object@tree$walks.force.labels)) {
tip.labels <- object@tree$walks.force.labels
} else {
# If there aren't any, generate some.
tip.labels <- treeForcePositionLabels(object)
}
# Choose short or long labels
if (use.short.names) tip.labels$label <- tip.labels$name.short else tip.labels$label <- tip.labels$name
# Limit to labels of cells that appear in the plot
segs.to.label <- intersect(unique(gg.data[which(gg.data$alpha > 0),"segment"]), tip.labels$seg)
tip.labels <- tip.labels[which(tip.labels$seg %in% segs.to.label),]
# Add the text
rgl::text3d(x=tip.labels$x, y=tip.labels$y, z=tip.labels$z, text=tip.labels$label, adj=0.5, cex=text.cex)
}
# Add a brief pause to ensure that rendering completes before moving on to next function.
Sys.sleep(0.1)
} else {
stop("Package rgl is required for this function. To install: install.packages('rgl')\n")
}
}
#' Determine position for segment labels in force-directed layout
#'
#' Finds cells near the end of each tip, and attempts to calculate a vector from
#' them to determine a good location for a label near each tip. Called automatically
#' by \code{\link{treeForceDirectedLayout}}, but can be re-run if needed.
#'
#' @param object An URD object
#' @param label.spacing (Numeric) Length of vector from final cell in each tip to label
#'
#' @return A data.frame with rows as tips, and columns containing labels and coordinates.
#' This is normally stored in \code{@@tree$walks.force.labels}.
#'
#' @examples
#' object@tree$walks.force.labels <- treeForcePositionLabels(object, label.spacing=2) # Move labels closer
#'
#' @export
treeForcePositionLabels <- function(object, label.spacing=5) {
# Grab the tips and their names
tip.labels <- data.frame(
seg=segTerminal(object),
name=object@tree$segment.names[segTerminal(object)],
name.short=object@tree$segment.names.short[segTerminal(object)],
stringsAsFactors=F, row.names=segTerminal(object)
)
# Anything that wasn't assigned a name previously: NA -> blank
tip.labels[is.na(tip.labels)] <- ""
# Get force-directed layout and add segment and pseudotime information
fdl.data <- object@tree$walks.force.layout
fdl.data$segment <- object@diff.data[rownames(fdl.data),"segment"]
fdl.data$pseudotime <- object@tree$pseudotime[rownames(fdl.data)]
# Order fdl data by pseudotime
fdl.data <- fdl.data[order(fdl.data$pseudotime, decreasing=T),]
# Figure out the location for each label by projecting out the vector at the end of each tip
tip.labels[,c("x","y","z")] <- t(as.data.frame(lapply(tip.labels$seg, function(tip) {
# Get cells from this tip
this.fdl <- fdl.data[which(fdl.data$segment==tip),]
# Find the vector describing direction that the cells at the tip are moving
tip.vec <- as.numeric(apply(this.fdl[c(1,5),c("x","y","telescope.pt")], 2, diff))
# Normalize the length of vector
s <- sqrt(label.spacing^2/(tip.vec[1]^2+tip.vec[2]^2+tip.vec[3]^2))
tip.vec <- s*tip.vec
# Place the label
label.coords <- apply(rbind(this.fdl[1,c("x","y","telescope.pt")], tip.vec),2,sum)
return(label.coords)
})))
return(tip.labels)
}
#' Calculate force-directed layout local density
#'
#' Add distance to nth nearest neighbor to force directed layout so that you can later
#' calculate local density and adjust alpha values.
#'
#' @importFrom RANN nn2
#'
#' @param object An URD object
#' @param neighbor (Numeric) Distance to \code{neighbor}th nearest neighbor is used for
#' density adjustment of point transparency on force-directed layout plots.
fdlDensity <- function(object, neighbor=10, knn=NULL) {
knn <- RANN::nn2(object@tree$walks.force.layout[,c("x","y","telescope.pt")], k = neighbor, treetype = 'kd')
object@tree$walks.force.layout$n.dist <- knn$nn.dists[,neighbor]
return(object)
}
#' Store a 3D view for force-directed layouts
#'
#' The last stored view will be the "default."
#'
#' @param object An URD object
#' @param view.name (Character) Name to use for this view
#'
#' @return An URD object with the parameters for that view stored in \code{@@tree$force.view.list}.
#'
#' @export
plotTreeForceStore3DView <- function(object, view.name) {
if (requireNamespace("rgl", quietly = TRUE)) {
# If not appending to a views list, make one
if (is.null(object@tree$force.view.list)) object@tree$force.view.list <- list()
# Grab rgl plot settings, if desired
rgl.settings <- rgl::par3d(no.readonly=F)
# Make a list of everything
this.view <- list(rgl.settings=rgl.settings)
# Store if in the views list
object@tree$force.view.list[[view.name]] <- this.view
# Store 'default' view
object@tree$force.view.default <- view.name
# Return the final views list
return(object)
} else {
stop("Package rgl is required for this function. To install: install.packages('rgl')\n")
}
}
#' Rotate coordinates in 3D
#'
#' @param m (Matrix) n cells x 3 matrix of coordinates
#' @param a (Numeric) Angle to rotate
#' @param axis (Character) Which axis to rotate around
#'
#' @return n x 3 matrix of rotated coordinates
#'
#' @keywords internal
rotateCoords3d <- function(m, a, axis=c("x","y","z")) {
if (length(axis) > 1) axis <- axis[1]
if (dim(m)[2] != 3) stop("Supply an n x 3 matrix.")
if (axis == "x") r.mat <- matrix(c(1, 0, 0, 0, cos(a), -1*sin(a), 0, sin(a), cos(a)), ncol=3, byrow=T)
if (axis == "y") r.mat <- matrix(c(cos(a), 0, sin(a), 0, 1, 0, -1*sin(a), 0, cos(a)), ncol=3, byrow=T)
if (axis == "z") r.mat <- matrix(c(cos(a), -1*sin(a), 0, sin(a), cos(a), 0, 0, 0, 1), ncol=3, byrow=T)
m %*% r.mat
}
#' Translate coordinates in 3D
#'
#' @param m (Matrix) n cells x 3 coordinates matrix
#' @param x (Numeric) Distance to move in x
#' @param y (Numeric) Distance to move in y
#' @param z (Numeric) Distance to move in z
#'
#' @return n x 3 matrix of translated coordinates
#'
#' @keywords internal
translateCoords3d <- function(m, x, y, z) {
if (dim(m)[2] != 3) stop("Supply an n x 3 matrix.")
sweep(m, 2, c(x, y, z), "+")
}
#' Rotate points in force-directed layout
#'
#' This function can be used in addition to \code{\link{treeForceTranslateCoords}}
#' in order to fine-tune the presentation of a force-directed layout. For instance,
#' this can be used to achieve an improved 2D visualization by modifying overlapping
#' portions of the layout.
#'
#' @param object An URD object
#' @param cells (Character vector) Cells to modify (Default \code{NULL} is all cells in the force-directed layout)
#' @param seg (Character) Instead of specifying cells, just grab all cells from this segment and downstream. Ignored if \code{cells} is specified.
#' @param angle (Numeric) Angle to rotate
#' @param axis (Character: "x", "y", "z") Axis around which to rotate
#' @param around.cell (Character or Numeric) If a character, then the name of a cell to rotate around. If numeric, then will choose the Nth cell by pseudotime. (This can be helpful for rotating a segment around its base to spread segments of the tree apart more. If \code{NULL}, then rotates around (0,0,0).
#' @param throw.out.cells (Numeric) Throw out the youngest N cells in pseudotime
#' @param pseudotime (Character) Pseudotime (i.e. column name of \code{@@pseudotime}) to use for determining cells to throw out.
#'
#' @return An URD object with the coordinates of some cells in \code{@@tree$walks.force.layout} modified.
#'
#' @export
treeForceRotateCoords <- function(object, cells=NULL, seg=NULL, angle, axis=c("x", "y", "z"), around.cell=NULL, throw.out.cells=0, pseudotime=NULL) {
# If segment specified, grab the cells from that segment
if (!is.null(seg) & is.null(cells)) cells <- cellsInCluster(object, "segment", segChildrenAll(object, seg, include.self=T))
# If cells not specified, do all cells in the layout.
if (is.null(cells)) cells <- rownames(object@tree$walks.force.layout)
# Make sure all cells specified are actually in the force-directed layout
cells <- intersect(cells, rownames(object@tree$walks.force.layout))
# If throw.out.cells is on, by removing n cells with the youngest pseudotime.
if ((throw.out.cells > 0)) {
if (is.null(pseudotime)) stop ("To use throw.out.cells, must specify pseudotime.")
cells <- setdiff(cells, cells[order(object@pseudotime[cells,pseudotime])[1:throw.out.cells]])
}
# Grab the appropriate chunk of the layout.
m <- as.matrix(object@tree$walks.force.layout[cells, c("x","y","telescope.pt")])
# If a cell is specified to rotate around, translate coords so that cell is at 0, 0, 0
if (!is.null(around.cell)) {
if (is.numeric(around.cell)) {
#around.cell <- rownames(m)[order(m[,3])[around.cell]]
around.cell <- cells[order(object@pseudotime[cells,pseudotime])[around.cell]]
}
around.coords <- m[around.cell,]
m <- translateCoords3d(m=m, x=-1*around.coords[1], y=-1*around.coords[2], z=-1*around.coords[3])
}
# Do the rotation
mprime <- rotateCoords3d(m=m, a=angle, axis=axis)
# If you were rotating around an arbitrary cell, put them back
if (!is.null(around.cell)) {
mprime <- translateCoords3d(m=mprime, x=around.coords[1], y=around.coords[2], z=around.coords[3])
}
object@tree$walks.force.layout[cells, c("x", "y", "telescope.pt")] <- mprime
return(object)
}
#' Translate points in force-directed layout
#'
#' This function can be used in addition to \code{\link{treeForceRotateCoords}}
#' in order to fine-tune the presentation of a force-directed layout. For instance,
#' this can be used to achieve an improved 2D visualization by modifying overlapping
#' portions of the layout.
#'
#' @param object An URD object
#' @param cells (Character vector) Cells to modify (Default \code{NULL} is all cells in the force-directed layout)
#' @param seg (Character) Instead of specifying cells, just grab all cells from this segment and downstream. Ignored if \code{cells} is specified.
#' @param x (Numeric) Distance to move cells along x-axis
#' @param y (Numeric) Distance to move cells along y-axis
#' @param z (Numeric) Distance to move cells along z-axis
#'
#' @return An URD object with the coordinates of some cells in \code{@@tree$walks.force.layout} modified.
#'
#' @export
treeForceTranslateCoords <- function(object, cells=NULL, seg=NULL, x, y, z) {
if (!is.null(seg) & is.null(cells)) cells <- cellsInCluster(object, "segment", segChildrenAll(object, seg, include.self=T))
if (is.null(cells)) cells <- rownames(object@tree$walks.force.layout)
cells <- intersect(cells, rownames(object@tree$walks.force.layout))
m <- as.matrix(object@tree$walks.force.layout[cells, c("x","y","telescope.pt")])
mprime <- translateCoords3d(m=m, x=x, y=y, z=z)
object@tree$walks.force.layout[cells, c("x", "y", "telescope.pt")] <- mprime
return(object)
}
#' Stretch points in force-directed layout
#'
#' This function can be used in addition to \code{\link{treeForceRotateCoords}}
#' and \code{\link{treeForceTranslateCoords}}
#' in order to fine-tune the presentation of a force-directed layout. For instance,
#' this can be used to achieve an improved 2D visualization by shortening or lengthening
#' a branch of the tree if branches have dramatically different lengths in pseudotime.
#' (This occurs, for instance, in data sets generated from a single timepoint in
#' a homeostatic animal instead of an explicit timecourse.) This transformation
#' will only occur in the force-directed layout, not the dendrogram tree layout.
#'
#' @param object An URD object
#' @param factor (Numeric) Factor to lengthen (>1) or shrink (<1) the branch by.
#' @param cells (Character vector) Cells to modify (Default \code{NULL} is all cells in the force-directed layout)
#' @param seg (Character) Instead of specifying cells, just grab all cells from this segment and downstream. Ignored if \code{cells} is specified.
#' @param throw.out.cells (Numeric) Throw out the youngest N cells in pseudotime (i.e. don't stretch those)
#'
#' @return An URD object with the coordinates of some cells in \code{@@tree$walks.force.layout} modified.
#'
#' @export
treeForceStretchCoords <- function(object, factor, cells=NULL, seg=NULL, throw.out.cells=0) {
if (!is.null(seg) & is.null(cells)) cells <- cellsInCluster(object, "segment", segChildrenAll(object, seg, include.self=T))
if (is.null(cells)) cells <- rownames(object@tree$walks.force.layout)
cells <- intersect(cells, rownames(object@tree$walks.force.layout))
if ((throw.out.cells > 0)) {
#if (is.null(pseudotime)) stop ("To use throw.out.cells, must specify pseudotime.")
cells <- setdiff(cells, cells[order(object@tree$pseudotime[cells])[1:throw.out.cells]])
#cells <- setdiff(cells, cells[order(object@pseudotime[cells,pseudotime])[1:throw.out.cells]])
}
m <- object@tree$walks.force.layout[cells, c("x","y","telescope.pt")]
m.min <- min(m$telescope.pt)
m$telescope.pt <- ((m$telescope.pt - m.min) * factor) + m.min
object@tree$walks.force.layout[cells, c("x", "y", "telescope.pt")] <- m
return(object)
}
treeForceStretchCoords <- function(object, factor, cells=NULL, seg=NULL, stretch.from=20, throw.out.cells=0, ramp=0.1) {
if (!is.null(seg) & is.null(cells)) cells <- cellsInCluster(object, "segment", segChildrenAll(object, seg, include.self=T))
if (is.null(cells)) cells <- rownames(object@tree$walks.force.layout)
cells <- intersect(cells, rownames(object@tree$walks.force.layout))
if ((throw.out.cells > 0)) {
#if (is.null(pseudotime)) stop ("To use throw.out.cells, must specify pseudotime.")
cells <- setdiff(cells, cells[order(object@tree$pseudotime[cells])[1:throw.out.cells]])
#cells <- setdiff(cells, cells[order(object@pseudotime[cells,pseudotime])[1:throw.out.cells]])
}
m <- object@tree$walks.force.layout[cells, c("x","y","telescope.pt")]
m.min <- sort(m$telescope.pt)[stretch.from]
m.max <- max(m$telescope.pt)
m.fac <- (m$telescope.pt - m.min) / (m.max - m.min)
m.fac <- pmax(m.fac, rep(0, length(m.fac)))
m.fac <- pmin(m.fac / ramp, rep(1, length(m.fac))) * factor
m$telescope.pt <- ((m$telescope.pt - m.min) * m.fac) + m.min
object@tree$walks.force.layout[cells, c("x", "y", "telescope.pt")] <- m
return(object)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.