From 229fea0ef25ac64ec9ddb9efcf6b3c216f5a9dd6 Mon Sep 17 00:00:00 2001 From: Michael Stadler Date: Sun, 21 Feb 2021 12:55:35 +0100 Subject: [PATCH] Add plotting functions to visualize RNA velocities (#30) * add plotVelocity * add return.intermediates argument to gridVectors * add tests for gridVectors(..., return.intermediates=TRUE) * add plotVelocityStream * fix minor issues and asthetics * add ggplot2 to Depends (required by metR::geom_streamline) * add metR to Imports * import all of ggplot2 * add tests for plotting functions * add tests for new utils * improve test coverage * move requirements for plotting to Suggests * do not plot or test without installed requirments * replace "colour" by "color" * use sym() in ggplot2::aes * move specific helpers to plotVelocity.R * remove or replace .isValidColor in plotVelocityStream * move .isValidColor from utils.R to plotVelocity.R * remove empty test file * remove or replace .isValidColor * remove redundant checks * use ggplot2/patchwork in plotVelocity instead of base graphics * add missing ggplot2:: * hide non-standard evaluation from R CMD check * dummy push to trigger new build on gha * use ggplot2::stat for double-dot aesthetics * dummy push to trigger new build on gha * add tests to improve coverage * fix NOTE about requireNamespace * use identical() instead of == when a single logical is expected Co-authored-by: Kevin Rue --- DESCRIPTION | 6 +- NAMESPACE | 5 + R/gridVectors.R | 23 ++- R/plotVelocity.R | 285 +++++++++++++++++++++++++++++++++ R/plotVelocityStream.R | 199 +++++++++++++++++++++++ man/gridVectors.Rd | 19 ++- man/plotVelocity.Rd | 105 ++++++++++++ man/plotVelocityStream.Rd | 121 ++++++++++++++ tests/testthat/test-embed.R | 4 + tests/testthat/test-plotting.R | 93 +++++++++++ 10 files changed, 855 insertions(+), 5 deletions(-) create mode 100644 R/plotVelocity.R create mode 100644 R/plotVelocityStream.R create mode 100644 man/plotVelocity.Rd create mode 100644 man/plotVelocityStream.Rd create mode 100644 tests/testthat/test-plotting.R diff --git a/DESCRIPTION b/DESCRIPTION index b1ec900..f4d76db 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -33,8 +33,12 @@ Suggests: scran, scater, scRNAseq, + Rtsne, + graphics, + grDevices, ggplot2, - Rtsne + patchwork, + metR StagedInstall: no License: MIT + file LICENSE Encoding: UTF-8 diff --git a/NAMESPACE b/NAMESPACE index d6bad73..a81e9a1 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -3,6 +3,8 @@ export(.grid_vectors) export(embedVelocity) export(gridVectors) +export(plotVelocity) +export(plotVelocityStream) export(scvelo) exportMethods(embedVelocity) exportMethods(gridVectors) @@ -22,6 +24,9 @@ importFrom(S4Vectors,selfmatch) importFrom(SingleCellExperiment,"reducedDim<-") importFrom(SingleCellExperiment,reducedDim) importFrom(SingleCellExperiment,reducedDimNames) +importFrom(SingleCellExperiment,reducedDims) +importFrom(SummarizedExperiment,assay) +importFrom(SummarizedExperiment,rowData) importFrom(basilisk,BasiliskEnvironment) importFrom(reticulate,import) importFrom(scuttle,normalizeCounts) diff --git a/R/gridVectors.R b/R/gridVectors.R index fccc711..65611d7 100644 --- a/R/gridVectors.R +++ b/R/gridVectors.R @@ -11,6 +11,8 @@ #' @param scale Logical scalar indicating whether the averaged vectors should be scaled by the grid resolution. #' @param as.data.frame Logical scalar indicating whether the output should be a data.frame. #' If \code{FALSE}, a list of two matrices is returned. +#' @param return.intermediates Logical scalar indicating whether intermediate objects should also be returned. +#' This enforces \code{as.data.frame=FALSE} and throws a warning if it is \code{TRUE}. #' @param ... For the generic, further arguments to pass to specific methods. #' #' For the SingleCellExperiment method, further arguments to pass to the ANY method. @@ -32,6 +34,13 @@ #' #' If \code{as.data.frame=TRUE}, a data.frame is returned with numeric columns of the same contents as the list above. #' Column names are prefixed by \code{start.*} and \code{end.*}. +#' +#' If \code{return.intermediates=TRUE}, a list is returned (irrespective of the value of \code{as.data.frame}) +#' that in addition to \code{start} and \code{end} also contains intermediate objects \code{limits} +#' (the ranges in x and y), \code{delta} (the grid intervals in x and y), \code{categories} (a +#' DataFrame with integer row and column indices for each cell that specify the grid field that it is +#' contained in), \code{grp} (numerical index of grid fields for each cell) and \code{vec} (velocity +#' vectors for non-empty grid fields). #' #' @author Aaron Lun #' @examples @@ -52,7 +61,7 @@ NULL #' @export #' @importFrom S4Vectors selfmatch DataFrame #' @importFrom stats median -.grid_vectors <- function(x, embedded, resolution=40, scale=TRUE, as.data.frame=TRUE) { +.grid_vectors <- function(x, embedded, resolution=40, scale=TRUE, as.data.frame=TRUE, return.intermediates=FALSE) { limits <- apply(x, 2, range) intercept <- limits[1,] delta <- (limits[2,] - intercept)/resolution @@ -72,8 +81,16 @@ NULL vec <- vec * target/median(d)/2 # divide by 2 to avoid running over into another block. } - FUN <- if (as.data.frame) data.frame else list - FUN(start=pos, end=pos + vec) + if (return.intermediates) { + if (as.data.frame) { + warning("ignoring 'as.data.frame=TRUE' for 'return.intermediates=TRUE'") + } + return(list(start=pos, end=pos + vec, limits=limits, delta=delta, + categories=categories, grp=grp, vec=vec)) + } else { + FUN <- if (as.data.frame) data.frame else list + return(FUN(start=pos, end=pos + vec)) + } } #' @export diff --git a/R/plotVelocity.R b/R/plotVelocity.R new file mode 100644 index 0000000..8d1ba24 --- /dev/null +++ b/R/plotVelocity.R @@ -0,0 +1,285 @@ +#' Phase and velocity graphs for a set of genes +#' +#' For a each gene in a set of genes, show the phase graph (spliced versus +#' unspliced counts and fitted model) and reduced dimension graphs with +#' cell colored by velocity and (spliced) expression. +#' +#' @param x A \linkS4class{SingleCellExperiment} object with RNA velocity results +#' as returned by \code{\link{scvelo}}, and low-dimensional coordinates, e.g., +#' after t-SNE, in its \code{\link{reducedDims}}. +#' @param genes A character vector with one or several genes for which to plot +#' phase and velocity graphs. \code{genes} have to be in \code{rownames(x)}. +#' @param use.dimred String or integer scalar specifying the reduced dimensions +#' to retrieve from \code{x}. +#' @param assay.splicedM An integer scalar or string specifying the assay of +#' \code{x} containing the moments of spliced abundances. +#' @param assay.unsplicedM An integer scalar or string specifying the assay of +#' \code{x} containing the moments unspliced abundances. +#' @param which.plots A character vector specifying which plots to create for +#' each gene. Possible values are \code{"phase", "velocity", "expression"} and +#' correspond to the phase graph or reduced dimension graphs with cells +#' colored by velocity or (spliced) expression. +#' @param genes.per.row An integer scalar with the numbers of genes to visualize +#' per row of plots. For example, if \code{which.plots = c("phase","expression")} +#' and \code{genes.per.row = 2}, the resulting figure will have four plot +#' panels per row. +#' @param color_by A character scalar specifying a column in \code{colData(x)} +#' to color cells in the phase graph. Alternatively, \code{color_by} can be +#' set to vector of valid R colors, either of length one (recycled for all +#' cells) or of length \code{ncol(x)}, which will then be used to color cells +#' in the phase graph. +#' @param color.alpha An integer scalar giving the transparency of colored +#' cells. Possible values are between 0 (fully transparent) and 1.0 (opaque). +#' @param colors.velocity,colors.expression Character vectors specifying the +#' color ranges used for mapping velocities and expression values. The +#' defaults are \code{RColorBrewer::brewer.pal(11, "RdYlGn")} for the +#' velocities and \code{viridisLite::viridis(11)} for the expression values. +#' @param max.abs.velo A numeric scalar greater than zero giving the maximum +#' absolute velocity to limit the color scale for the \code{"velocity"} graph. +#' +#' @details Please note that \code{plotVelocity} will modify parameters of +#' the current graphics device using \code{\link{layout}} and \code{\link{par}}, +#' in order to create the layout for the generated graph panels. +#' +#' @return A patchwork object with the plots selected by \code{which.plot} for +#' the genes in \code{genes}, arranged in a grid according to \code{genes.per.row}. +#' +#' @author Michael Stadler +#' +#' @examples +#' library(scuttle) +#' set.seed(42) +#' sce1 <- mockSCE(ncells = 100, ngenes = 500) +#' sce2 <- mockSCE(ncells = 100, ngenes = 500) +#' +#' datlist <- list(X=counts(sce1), spliced=counts(sce1), unspliced=counts(sce2)) +#' +#' out1 <- scvelo(datlist, mode = "steady_state") +#' out2 <- scvelo(datlist, mode = "dynamical") +#' +#' plotVelocity(out1, c("Gene_0031","Gene_0268")) +#' plotVelocity(out2, c("Gene_0031","Gene_0268")) +#' +#' @seealso +#' \code{\link{scvelo}}, to generate \code{x}, +#' \code{\link[RColorBrewer]{brewer.pal}} and \code{\link[viridisLite]{viridis}} +#' for creation of color palettes, packages \pkg{ggplot2} and \pkg{patchwork} +#' used to generate and arrange the plots. +#' +#' @export +#' @importFrom SummarizedExperiment assay rowData +#' @importFrom SingleCellExperiment reducedDim reducedDims reducedDimNames +plotVelocity <- function(x, genes, use.dimred = 1, + assay.splicedM = "Ms", assay.unsplicedM = "Mu", + which.plots = c("phase", "velocity", "expression"), + genes.per.row = 1, + color_by = "#222222", + color.alpha = 0.4, + colors.velocity = c("#A50026", "#D73027", "#F46D43", + "#FDAE61", "#FEE08B", "#FFFFBF", + "#D9EF8B", "#A6D96A", "#66BD63", + "#1A9850", "#006837"), + colors.expression = c("#440154", "#482576", "#414487", + "#35608D", "#2A788E", "#21908C", + "#22A884", "#43BF71", "#7AD151", + "#BBDF27", "#FDE725"), + max.abs.velo = 0.001) { + # check arguments + genes <- unique(genes) + if (!all(genes %in% rownames(x))) { + stop("not all 'genes' are not found in 'x'") + } + stopifnot(exprs = { + all(which.plots %in% c("phase", "velocity", "expression")) + is.numeric(genes.per.row) + length(genes.per.row) == 1L + is.numeric(max.abs.velo) + length(max.abs.velo) == 1L + max.abs.velo >= 0.0 + }) + if (is.numeric(use.dimred)) { + stopifnot(exprs = { + length(use.dimred) == 1L + use.dimred <= length(reducedDims(x)) + }) + use.dimred <- reducedDimNames(x)[use.dimred] + } else if (is.character(use.dimred)) { + stopifnot(exprs = { + length(use.dimred) == 1L + use.dimred %in% reducedDimNames(x) + }) + } else { + stop("'use.dimred' is not a valid value for use in reducedDim(x, use.dimred)") + } + + # extract data from x + S <- assay(x, assay.splicedM) + U <- assay(x, assay.unsplicedM) + V <- assay(x, "velocity") + rd <- rowData(x) + xy <- reducedDim(x, use.dimred) + df2 <- data.frame(x = xy[, 1], + y = xy[, 2]) + + + # iterate over genes and create graphs + pL <- list() + for (gene in genes) { + df1 <- data.frame(s = as.vector(S[gene, ]), + u = as.vector(U[gene, ])) + + # spliced/unspliced phase portrait with model estimates + if ("phase" %in% which.plots) { + if (is.character(color_by) && length(color_by) == 1L && color_by %in% colnames(colData(x))) { + df1$col <- factor(colData(x)[, color_by]) + colByFeat <- TRUE + } else if (length(color_by) == 1L || length(color_by) == ncol(x)) { + df1$col <- color_by + colByFeat <- FALSE + } else { + stop("invalid 'colour_by' (not in colData(x), nor of the right length)") + } + p1 <- ggplot2::ggplot(df1, ggplot2::aes(x = !!ggplot2::sym("s"), + y = !!ggplot2::sym("u"))) + + ggplot2::xlim(0, max(df1$s, na.rm = TRUE)) + + ggplot2::ylim(0, max(df1$u, na.rm = TRUE)) + + ggplot2::labs(title = gene, x = "spliced", y = "unspliced") + + ggplot2::theme_minimal() + + ggplot2::theme(panel.grid.major = ggplot2::element_blank(), + panel.grid.minor = ggplot2::element_blank(), + axis.ticks = ggplot2::element_line(color = "grey20"), + panel.border = ggplot2::element_rect(fill = NA, color = "grey20")) + if (colByFeat) { + p1 <- p1 + ggplot2::geom_point(ggplot2::aes(color = !!ggplot2::sym("col")), + alpha = color.alpha, show.legend = FALSE) + } else { + p1 <- p1 + ggplot2::geom_point(color = df1$col, alpha = color.alpha) + } + + # gene specific model parameters + # scvelo(..., mode = "dynamical") creates "fit_" prefixes, + # all other modes only have "velocity_gamma" + if ("fit_gamma" %in% colnames(rd)) { + alpha <- rd[gene, "fit_alpha"] + beta <- rd[gene, "fit_beta"] * rd[gene, "fit_scaling"] + gamma <- rd[gene, "fit_gamma"] + scaling <- rd[gene, "fit_scaling"] + # t_ is the pseudotime separating the increase/decrease phases + t_ <- rd[gene, "fit_t_"] + } else { + alpha <- beta <- scaling <- t_ <- NULL + if ("velocity_gamma" %in% colnames(rd)) { + gamma <- rd[gene, "velocity_gamma"] + } else { + gamma <- NULL + } + } + + # steady-state ratio + if (!is.null(gamma) && is.finite(gamma)) { + if (!is.null(beta) && is.finite(beta)) { + # in the dynamical model + p1 <- p1 + ggplot2::geom_abline(intercept = 0, + slope = gamma / beta * scaling, + linetype = "dashed") + } else { + # in a static model + p1 <- p1 + ggplot2::geom_abline(intercept = 0, + slope = gamma, + linetype = "dashed") + } + } + + # dynamic model curve + if (!is.null(beta) && is.finite(beta)) { + # calculate spliced and unspliced abundances from initial + # conditions and fitted model parameters + # see https://github.com/theislab/scvelo/blob/95d90de3d0935ce58a01218c9f179c9494ff593e/scvelo/plotting/simulation.py#L26 + if (all(c("fit_u0","fit_s0") %in% colnames(rd))) { + u0_offset <- rd[gene, "fit_u0"] + s0_offset <- rd[gene, "fit_s0"] + } else { + u0_offset <- s0_offset <- 0 + } + # see https://github.com/theislab/scvelo/blob/95d90de3d0935ce58a01218c9f179c9494ff593e/scvelo/tools/dynamical_model_utils.py#L602 + ft <- as.vector(assay(x, "fit_t")[gene, ]) + o <- as.numeric(ft < t_) # 1 or 0 for cells increasing or decreasing this gene's expression + tau <- ft * o + (ft - t_) * (1 - o) + u0 <- s0 <- alpha0 <- 0 + u0_ <- .unspliced(t_, u0, alpha, beta) + s0_ <- .spliced(t_, s0, u0, alpha, beta, gamma) + u0 <- u0 * o + u0_ * (1 - o) + s0 <- s0 * o + s0_ * (1 - o) + alpha <- alpha * o + alpha0 * (1 - o) + # see https://github.com/theislab/scvelo/blob/95d90de3d0935ce58a01218c9f179c9494ff593e/scvelo/plotting/simulation.py#L54 + ut <- .unspliced(tau, u0, alpha, beta) * scaling + u0_offset + st <- .spliced(tau, s0, u0, alpha, beta, gamma) + s0_offset + p1 <- p1 + ggplot2::geom_point(data = data.frame(x = st, y = ut), + mapping = ggplot2::aes(x = !!ggplot2::sym("x"), + y = !!ggplot2::sym("y")), + color = "purple", shape = 46) + } + pL <- c(pL, list(p1)) + } + + # velocity plot + if ("velocity" %in% which.plots) { + df2$z <- as.vector(V[gene, ]) + if (all(!is.finite(df2$z))) { + warning("velocity estimates for ", gene, " are all non-finite") + } + mx <- max(c(max.abs.velo, abs(df2$z)), na.rm = TRUE) + p2 <- ggplot2::ggplot(df2, ggplot2::aes(x = !!ggplot2::sym("x"), + y = !!ggplot2::sym("y"), + color = !!ggplot2::sym("z"))) + + ggplot2::geom_point(alpha = color.alpha) + + ggplot2::scale_color_gradientn(colours = colors.velocity, + limits = c(-mx, mx)) + + ggplot2::labs(title = "velocity", x = paste0(use.dimred, " 1"), + y = paste0(use.dimred, " 2"), color = ggplot2::element_blank()) + + ggplot2::theme_minimal() + + ggplot2::theme(axis.text = ggplot2::element_blank(), + panel.grid.major = ggplot2::element_blank(), + panel.grid.minor = ggplot2::element_blank()) + pL <- c(pL, list(p2)) + } + + # expression plot + if ("expression" %in% which.plots) { + df2$z <- as.vector(S[gene, ]) + p3 <- ggplot2::ggplot(df2, ggplot2::aes(x = !!ggplot2::sym("x"), + y = !!ggplot2::sym("y"), + color = !!ggplot2::sym("z"))) + + ggplot2::geom_point(alpha = color.alpha) + + ggplot2::scale_color_gradientn(colours = colors.expression) + + ggplot2::labs(title = "expression", x = paste0(use.dimred, " 1"), + y = paste0(use.dimred, " 2"), color = ggplot2::element_blank()) + + ggplot2::theme_minimal() + + ggplot2::theme(axis.text = ggplot2::element_blank(), + panel.grid.major = ggplot2::element_blank(), + panel.grid.minor = ggplot2::element_blank()) + pL <- c(pL, list(p3)) + } + } + + # create plot panel layout + nplts <- length(unique(which.plots)) + nrows <- ceiling(length(genes) / genes.per.row) + p <- do.call(patchwork::wrap_plots, c(pL, list(nrow = nrows))) + + return(p) +} + +# helper functions to calculate (un)spliced counts from starting point and model parameters +# see https://github.com/theislab/scvelo/blob/95d90de3d0935ce58a01218c9f179c9494ff593e/scvelo/tools/dynamical_model_utils.py#L109 +.unspliced <- function(tau, u0, alpha, beta) { + expu <- exp(-beta * tau) + return(u0 * expu + alpha / beta * (1 - expu)) +} +# see https://github.com/theislab/scvelo/blob/95d90de3d0935ce58a01218c9f179c9494ff593e/scvelo/tools/dynamical_model_utils.py#L114 +.spliced <- function(tau, s0, u0, alpha, beta, gamma) { + cn <- (alpha - u0 * beta) * 1 / (gamma - beta) + expu <- exp(-beta * tau) + exps <- exp(-gamma * tau) + return(s0 * exps + alpha / gamma * (1 - exps) + cn * (exps - expu)) +} diff --git a/R/plotVelocityStream.R b/R/plotVelocityStream.R new file mode 100644 index 0000000..fcb473e --- /dev/null +++ b/R/plotVelocityStream.R @@ -0,0 +1,199 @@ +#' Velocity stream plot in low-dimensional space +#' +#' Plot velocities embedded into low-dimensional space as a stream plot. Stream +#' lines are lines that follow the gradient in the velocity field and illustrate +#' paths that cells could follow based on observed RNA velocities. +#' +#' @param sce A \linkS4class{SingleCellExperiment} object containing +#' low-dimensional coordinates, e.g., after t-SNE, in its +#' \code{\link{reducedDims}}. +#' @param embedded A low-dimensional projection of the velocity vectors into the +#' embedding of \code{sce}. This should be of the same dimensions as \code{sce} +#' and is typically produced by \code{\link{embedVelocity}}. +#' @param use.dimred String or integer scalar specifying the reduced dimensions +#' to retrieve from \code{sce}. +#' @param color_by A character scalar specifying a column in \code{colData(sce)} +#' to color cells in the phase graph. Alternatively, \code{color_by} can be +#' set to a valid R color to be used to color cells. +#' @param color.alpha An integer scalar giving the transparency of colored +#' cells. Possible values are between 0 (fully transparent) and 1.0 (opaque). +#' @param grid.resolution Integer scalar specifying the resolution of the grid, +#' in terms of the number of grid intervals along each axis. +#' @param scale Logical scalar indicating whether the averaged vectors should be +#' scaled by the grid resolution. +#' @param stream.L Integer scalar giving the typical length of a streamline +#' low-dimensional space units. +#' @param stream.min.L A numeric scalar with the minimum length of segments to be shown. +#' @param stream.res Numeric scalar specifying the resolution of estimated +#' streamlines (higher numbers increase smoothness of lines but also the time +#' for computation). +#' @param stream.width A numeric scalar controlling the width of streamlines. +#' @param color.streamlines Logical scalar. If \code{TRUE} streamlines will +#' be colored by local velocity. Arrows cannot be shown in that case. +#' @param color.streamlines.map A character vector specifying the +#' color range used for mapping local velocities to streamline colors. The +#' default is \code{viridisLite::viridis(11)}. +#' @param arrow.angle,arrow.length Numeric scalars giving the \code{angle} and +#' \code{length} of arrowheads. +#' +#' @details \code{grid.resolution} and \code{scale} are passed to +#' \code{\link{gridVectors}}, which is used to summarized the velocity vectors +#' into an initial grid. A full regular grid is computed from that and used +#' in \code{\link[metR]{geom_streamline}} to calculate streamlines. The +#' following arguments are passed to the arguments given in parenthesis of +#' \code{\link[metR]{geom_streamline}}: +#' \code{stream.L} (\code{L}), \code{stream.res} (\code{res}), +#' \code{stream.min.L} (\code{min.L}), \code{arrow.angle} (\code{arrow.angle}) +#' and \code{arrow.length} (\code{arrow.length}). +#' Streamlines are computed by simple integration with a forward Euler method, +#' and \code{stream.L} and \code{stream.res} are used to compute the number of +#' steps and the time interval between steps for the integration. +#' \code{stream.width} is multiplied with \code{..step..} estimated by +#' \code{\link[metR]{geom_streamline}} to control the width of streamlines. +#' +#' @return A \code{ggplot2} object with the streamline plot. +#' +#' @author Michael Stadler +#' +#' @examples +#' library(scuttle) +#' set.seed(42) +#' sce1 <- mockSCE(ncells = 100, ngenes = 500) +#' sce2 <- mockSCE(ncells = 100, ngenes = 500) +#' +#' datlist <- list(X=counts(sce1), spliced=counts(sce1), unspliced=counts(sce2)) +#' +#' out <- scvelo(datlist, mode = "dynamical") +#' +#' em <- embedVelocity(reducedDim(out, 1), out)[,1:2] +#' +#' plotVelocityStream(out, em) +#' plotVelocityStream(out, em, color.streamlines = TRUE) +#' +#' @seealso \code{\link{gridVectors}} used to summarize velocity vectors into +#' a grid (velocity field), the \pkg{ggplot2} package used for plotting, +#' \code{\link[metR]{geom_streamline}} in package \pkg{metR} used to +#' calculate and add streamlines from the RNA velocity field to the plot, +#' \code{\link[viridisLite]{viridis}} for creation of color palettes. +#' +#' @export +#' @importFrom S4Vectors DataFrame +plotVelocityStream <- function(sce, embedded, use.dimred = 1, + color_by = "#444444", color.alpha = 0.2, + grid.resolution = 60, scale = TRUE, + stream.L = 10, stream.min.L = 0, stream.res = 4, + stream.width = 8, + color.streamlines = FALSE, + color.streamlines.map = c("#440154", "#482576", "#414487", + "#35608D", "#2A788E", "#21908C", + "#22A884", "#43BF71", "#7AD151", + "#BBDF27", "#FDE725"), + arrow.angle = 8, arrow.length = 0.8) { + if (!identical(ncol(sce), nrow(embedded))) { + stop("'sce' and 'embedded' do not have consistent dimensions.") + } + if (is.numeric(use.dimred)) { + stopifnot(exprs = { + identical(length(use.dimred), 1L) + use.dimred <= length(reducedDims(sce)) + }) + use.dimred <- reducedDimNames(sce)[use.dimred] + } + else if (is.character(use.dimred)) { + stopifnot(exprs = { + length(use.dimred) == 1L + use.dimred %in% reducedDimNames(sce) + }) + } + else { + stop("'use.dimred' is not a valid value for use in reducedDim(sce, use.dimred)") + } + if (!requireNamespace("ggplot2")) { + stop("'plotVelocityStream' requires the package 'ggplot2'.") + } + + # get coordinates in reduced dimensional space + xy <- reducedDim(sce, use.dimred)[, 1:2] + + # summarize velocities in a grid + gr <- gridVectors(x = xy, embedded = embedded, + resolution = grid.resolution, scale = scale, + as.data.frame = FALSE, + return.intermediates = TRUE) + + # now make it a regular grid needed for metR::geom_streamline + xbreaks <- seq(gr$limits[1,1], gr$limits[2,1], by = gr$delta[1]) + ybreaks <- seq(gr$limits[1,2], gr$limits[2,2], by = gr$delta[2]) + plotdat2 <- expand.grid(x = xbreaks + gr$delta[1] / 2, + y = ybreaks + gr$delta[2] / 2, + dx = 0, dy = 0) + allcategories <- DataFrame(expand.grid(V1 = seq(0, grid.resolution), + V2 = seq(0, grid.resolution))) + ivec <- match(gr$categories[sort(unique(gr$grp)), ], allcategories) + plotdat2[ivec, c("dx", "dy")] <- gr$vec + + + # plot it using ggplot2 and metR::geom_streamline + plotdat1 <- data.frame(xy) + colnames(plotdat1) <- c("x", "y") + if (is.character(color_by) && length(color_by) == 1L && color_by %in% colnames(colData(sce))) { + plotdat1 <- cbind(plotdat1, col = colData(sce)[, color_by]) + colByFeat <- TRUE + } else { + colByFeat <- FALSE + } + p <- ggplot2::ggplot(plotdat1, ggplot2::aes(x = !!ggplot2::sym("x"), y = !!ggplot2::sym("y"))) + + ggplot2::labs(x = paste(use.dimred, "1"), y = paste(use.dimred, "2")) + if (!colByFeat) { + colMatrix <- grDevices::col2rgb(col = color_by, alpha = TRUE) + if (any(colMatrix[4, ] != 255)) { + warning("ignoring 'color.alpha' as 'color_by' already specifies alpha channels") + color.alpha <- colMatrix[4, ] / 255 + } + p <- p + ggplot2::geom_point(color = color_by, alpha = color.alpha) + } else { + p <- p + ggplot2::geom_point(ggplot2::aes(color = !!ggplot2::sym("col")), alpha = color.alpha) + + ggplot2::labs(color = color_by) + } + if (color.streamlines) { + # remark: when coloring streamlines, we currently cannot have any arrows + # remark: ..dx.., ..dy.. and ..step.. are calculated by metR::geom_streamline + p <- p + + metR::geom_streamline(mapping = ggplot2::aes(x = !!ggplot2::sym("x"), + y = !!ggplot2::sym("y"), + dx = !!ggplot2::sym("dx"), + dy = !!ggplot2::sym("dy"), + size = stream.width * !!ggplot2::sym("..step.."), + alpha = !!ggplot2::sym("..step.."), + color = ggplot2::stat(sqrt((!!ggplot2::sym("..dx.."))^2 + + (!!ggplot2::sym("..dy.."))^2))), + arrow = NULL, lineend = "round", + data = plotdat2, size = 0.6, jitter = 2, + L = stream.L, min.L = stream.min.L, + res = stream.res, inherit.aes = FALSE) + + ggplot2::scale_color_gradientn(colors = color.streamlines.map, + guide = "none") + + ggplot2::scale_alpha_continuous(guide = "none") + + ggplot2::theme_minimal() + + ggplot2::theme(axis.text = ggplot2::element_blank(), + panel.grid.major = ggplot2::element_blank(), + panel.grid.minor = ggplot2::element_blank()) + } else { + p <- p + + metR::geom_streamline(mapping = ggplot2::aes(x = !!ggplot2::sym("x"), + y = !!ggplot2::sym("y"), + dx = !!ggplot2::sym("dx"), + dy = !!ggplot2::sym("dy"), + size = stream.width * !!ggplot2::sym("..step..")), + data = plotdat2, size = 0.3, jitter = 2, + L = stream.L, min.L = stream.min.L, + res = stream.res, arrow.angle = arrow.angle, + arrow.length = arrow.length, inherit.aes = FALSE) + + ggplot2::theme_minimal() + + ggplot2::theme(axis.text = ggplot2::element_blank(), + panel.grid.major = ggplot2::element_blank(), + panel.grid.minor = ggplot2::element_blank()) + } + + return(p) +} diff --git a/man/gridVectors.Rd b/man/gridVectors.Rd index 0502d2d..777eaef 100644 --- a/man/gridVectors.Rd +++ b/man/gridVectors.Rd @@ -8,7 +8,14 @@ \usage{ gridVectors(x, embedded, ...) -\S4method{gridVectors}{ANY}(x, embedded, resolution = 40, scale = TRUE, as.data.frame = TRUE) +\S4method{gridVectors}{ANY}( + x, + embedded, + resolution = 40, + scale = TRUE, + as.data.frame = TRUE, + return.intermediates = FALSE +) \S4method{gridVectors}{SingleCellExperiment}(x, embedded, ..., use.dimred = 1) } @@ -31,6 +38,9 @@ in terms of the number of grid intervals along each axis.} \item{as.data.frame}{Logical scalar indicating whether the output should be a data.frame. If \code{FALSE}, a list of two matrices is returned.} +\item{return.intermediates}{Logical scalar indicating whether intermediate objects should also be returned. +This enforces \code{as.data.frame=FALSE} and throws a warning if it is \code{TRUE}.} + \item{use.dimred}{String or integer scalar specifying the reduced dimensions to retrieve from \code{x}.} } \value{ @@ -41,6 +51,13 @@ and \code{end} contains the endpoint after adding the (scaled) average of the bl If \code{as.data.frame=TRUE}, a data.frame is returned with numeric columns of the same contents as the list above. Column names are prefixed by \code{start.*} and \code{end.*}. + +If \code{return.intermediates=TRUE}, a list is returned (irrespective of the value of \code{as.data.frame}) +that in addition to \code{start} and \code{end} also contains intermediate objects \code{limits} +(the ranges in x and y), \code{delta} (the grid intervals in x and y), \code{categories} (a +DataFrame with integer row and column indices for each cell that specify the grid field that it is +contained in), \code{grp} (numerical index of grid fields for each cell) and \code{vec} (velocity +vectors for non-empty grid fields). } \description{ Summarize the velocity vectors into a grid, usually for easy plotting. diff --git a/man/plotVelocity.Rd b/man/plotVelocity.Rd new file mode 100644 index 0000000..0a04d10 --- /dev/null +++ b/man/plotVelocity.Rd @@ -0,0 +1,105 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plotVelocity.R +\name{plotVelocity} +\alias{plotVelocity} +\title{Phase and velocity graphs for a set of genes} +\usage{ +plotVelocity( + x, + genes, + use.dimred = 1, + assay.splicedM = "Ms", + assay.unsplicedM = "Mu", + which.plots = c("phase", "velocity", "expression"), + genes.per.row = 1, + color_by = "#222222", + color.alpha = 0.4, + colors.velocity = c("#A50026", "#D73027", "#F46D43", "#FDAE61", "#FEE08B", "#FFFFBF", + "#D9EF8B", "#A6D96A", "#66BD63", "#1A9850", "#006837"), + colors.expression = c("#440154", "#482576", "#414487", "#35608D", "#2A788E", + "#21908C", "#22A884", "#43BF71", "#7AD151", "#BBDF27", "#FDE725"), + max.abs.velo = 0.001 +) +} +\arguments{ +\item{x}{A \linkS4class{SingleCellExperiment} object with RNA velocity results +as returned by \code{\link{scvelo}}, and low-dimensional coordinates, e.g., +after t-SNE, in its \code{\link{reducedDims}}.} + +\item{genes}{A character vector with one or several genes for which to plot +phase and velocity graphs. \code{genes} have to be in \code{rownames(x)}.} + +\item{use.dimred}{String or integer scalar specifying the reduced dimensions +to retrieve from \code{x}.} + +\item{assay.splicedM}{An integer scalar or string specifying the assay of +\code{x} containing the moments of spliced abundances.} + +\item{assay.unsplicedM}{An integer scalar or string specifying the assay of +\code{x} containing the moments unspliced abundances.} + +\item{which.plots}{A character vector specifying which plots to create for +each gene. Possible values are \code{"phase", "velocity", "expression"} and +correspond to the phase graph or reduced dimension graphs with cells +colored by velocity or (spliced) expression.} + +\item{genes.per.row}{An integer scalar with the numbers of genes to visualize +per row of plots. For example, if \code{which.plots = c("phase","expression")} +and \code{genes.per.row = 2}, the resulting figure will have four plot +panels per row.} + +\item{color_by}{A character scalar specifying a column in \code{colData(x)} +to color cells in the phase graph. Alternatively, \code{color_by} can be +set to vector of valid R colors, either of length one (recycled for all +cells) or of length \code{ncol(x)}, which will then be used to color cells +in the phase graph.} + +\item{color.alpha}{An integer scalar giving the transparency of colored +cells. Possible values are between 0 (fully transparent) and 1.0 (opaque).} + +\item{colors.velocity, colors.expression}{Character vectors specifying the +color ranges used for mapping velocities and expression values. The +defaults are \code{RColorBrewer::brewer.pal(11, "RdYlGn")} for the +velocities and \code{viridisLite::viridis(11)} for the expression values.} + +\item{max.abs.velo}{A numeric scalar greater than zero giving the maximum +absolute velocity to limit the color scale for the \code{"velocity"} graph.} +} +\value{ +A patchwork object with the plots selected by \code{which.plot} for + the genes in \code{genes}, arranged in a grid according to \code{genes.per.row}. +} +\description{ +For a each gene in a set of genes, show the phase graph (spliced versus +unspliced counts and fitted model) and reduced dimension graphs with +cell colored by velocity and (spliced) expression. +} +\details{ +Please note that \code{plotVelocity} will modify parameters of + the current graphics device using \code{\link{layout}} and \code{\link{par}}, + in order to create the layout for the generated graph panels. +} +\examples{ +library(scuttle) +set.seed(42) +sce1 <- mockSCE(ncells = 100, ngenes = 500) +sce2 <- mockSCE(ncells = 100, ngenes = 500) + +datlist <- list(X=counts(sce1), spliced=counts(sce1), unspliced=counts(sce2)) + +out1 <- scvelo(datlist, mode = "steady_state") +out2 <- scvelo(datlist, mode = "dynamical") + +plotVelocity(out1, c("Gene_0031","Gene_0268")) +plotVelocity(out2, c("Gene_0031","Gene_0268")) + +} +\seealso{ +\code{\link{scvelo}}, to generate \code{x}, +\code{\link[RColorBrewer]{brewer.pal}} and \code{\link[viridisLite]{viridis}} +for creation of color palettes, packages \pkg{ggplot2} and \pkg{patchwork} +used to generate and arrange the plots. +} +\author{ +Michael Stadler +} diff --git a/man/plotVelocityStream.Rd b/man/plotVelocityStream.Rd new file mode 100644 index 0000000..c9fa988 --- /dev/null +++ b/man/plotVelocityStream.Rd @@ -0,0 +1,121 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plotVelocityStream.R +\name{plotVelocityStream} +\alias{plotVelocityStream} +\title{Velocity stream plot in low-dimensional space} +\usage{ +plotVelocityStream( + sce, + embedded, + use.dimred = 1, + color_by = "#444444", + color.alpha = 0.2, + grid.resolution = 60, + scale = TRUE, + stream.L = 10, + stream.min.L = 0, + stream.res = 4, + stream.width = 8, + color.streamlines = FALSE, + color.streamlines.map = c("#440154", "#482576", "#414487", "#35608D", "#2A788E", + "#21908C", "#22A884", "#43BF71", "#7AD151", "#BBDF27", "#FDE725"), + arrow.angle = 8, + arrow.length = 0.8 +) +} +\arguments{ +\item{sce}{A \linkS4class{SingleCellExperiment} object containing +low-dimensional coordinates, e.g., after t-SNE, in its +\code{\link{reducedDims}}.} + +\item{embedded}{A low-dimensional projection of the velocity vectors into the +embedding of \code{sce}. This should be of the same dimensions as \code{sce} +and is typically produced by \code{\link{embedVelocity}}.} + +\item{use.dimred}{String or integer scalar specifying the reduced dimensions +to retrieve from \code{sce}.} + +\item{color_by}{A character scalar specifying a column in \code{colData(sce)} +to color cells in the phase graph. Alternatively, \code{color_by} can be +set to a valid R color to be used to color cells.} + +\item{color.alpha}{An integer scalar giving the transparency of colored +cells. Possible values are between 0 (fully transparent) and 1.0 (opaque).} + +\item{grid.resolution}{Integer scalar specifying the resolution of the grid, +in terms of the number of grid intervals along each axis.} + +\item{scale}{Logical scalar indicating whether the averaged vectors should be +scaled by the grid resolution.} + +\item{stream.L}{Integer scalar giving the typical length of a streamline +low-dimensional space units.} + +\item{stream.min.L}{A numeric scalar with the minimum length of segments to be shown.} + +\item{stream.res}{Numeric scalar specifying the resolution of estimated +streamlines (higher numbers increase smoothness of lines but also the time +for computation).} + +\item{stream.width}{A numeric scalar controlling the width of streamlines.} + +\item{color.streamlines}{Logical scalar. If \code{TRUE} streamlines will +be colored by local velocity. Arrows cannot be shown in that case.} + +\item{color.streamlines.map}{A character vector specifying the +color range used for mapping local velocities to streamline colors. The +default is \code{viridisLite::viridis(11)}.} + +\item{arrow.angle, arrow.length}{Numeric scalars giving the \code{angle} and +\code{length} of arrowheads.} +} +\value{ +A \code{ggplot2} object with the streamline plot. +} +\description{ +Plot velocities embedded into low-dimensional space as a stream plot. Stream +lines are lines that follow the gradient in the velocity field and illustrate +paths that cells could follow based on observed RNA velocities. +} +\details{ +\code{grid.resolution} and \code{scale} are passed to + \code{\link{gridVectors}}, which is used to summarized the velocity vectors + into an initial grid. A full regular grid is computed from that and used + in \code{\link[metR]{geom_streamline}} to calculate streamlines. The + following arguments are passed to the arguments given in parenthesis of + \code{\link[metR]{geom_streamline}}: + \code{stream.L} (\code{L}), \code{stream.res} (\code{res}), + \code{stream.min.L} (\code{min.L}), \code{arrow.angle} (\code{arrow.angle}) + and \code{arrow.length} (\code{arrow.length}). + Streamlines are computed by simple integration with a forward Euler method, + and \code{stream.L} and \code{stream.res} are used to compute the number of + steps and the time interval between steps for the integration. + \code{stream.width} is multiplied with \code{..step..} estimated by + \code{\link[metR]{geom_streamline}} to control the width of streamlines. +} +\examples{ +library(scuttle) +set.seed(42) +sce1 <- mockSCE(ncells = 100, ngenes = 500) +sce2 <- mockSCE(ncells = 100, ngenes = 500) + +datlist <- list(X=counts(sce1), spliced=counts(sce1), unspliced=counts(sce2)) + +out <- scvelo(datlist, mode = "dynamical") + +em <- embedVelocity(reducedDim(out, 1), out)[,1:2] + +plotVelocityStream(out, em) +plotVelocityStream(out, em, color.streamlines = TRUE) + +} +\seealso{ +\code{\link{gridVectors}} used to summarize velocity vectors into + a grid (velocity field), the \pkg{ggplot2} package used for plotting, + \code{\link[metR]{geom_streamline}} in package \pkg{metR} used to + calculate and add streamlines from the RNA velocity field to the plot, + \code{\link[viridisLite]{viridis}} for creation of color palettes. +} +\author{ +Michael Stadler +} diff --git a/tests/testthat/test-embed.R b/tests/testthat/test-embed.R index 0c4e0c9..e1cadd0 100644 --- a/tests/testthat/test-embed.R +++ b/tests/testthat/test-embed.R @@ -29,6 +29,10 @@ test_that("gridVectors works correctly", { out <- gridVectors(tsne.results, projected) expect_identical(ncol(out), 4L) + expect_warning(out <- gridVectors(tsne.results, projected, return.intermediates=TRUE)) + expect_type(out, "list") + expect_length(out, 7L) + out <- gridVectors(tsne.results, projected, as.data.frame=FALSE) expect_type(out, "list") diff --git a/tests/testthat/test-plotting.R b/tests/testthat/test-plotting.R new file mode 100644 index 0000000..c829932 --- /dev/null +++ b/tests/testthat/test-plotting.R @@ -0,0 +1,93 @@ +# This tests the plotting functionality. +# library(testthat); library(velociraptor); source("test-plotting.R") + +library(scuttle) +#' set.seed(42) +sce1 <- mockSCE() +sce2 <- mockSCE() + +spliced <- counts(sce1) +unspliced <- counts(sce2) +datlist <- list(X=counts(sce1), spliced=counts(sce1), unspliced=counts(sce2)) + +set.seed(100001) +out1 <- scvelo(datlist, mode = "steady_state") +out2 <- scvelo(datlist, mode = "dynamical") + +genes1 <- rownames(out2)[!is.finite(rowData(out2)$fit_alpha)][1:2] +genes2 <- rownames(out2)[is.finite(rowData(out2)$fit_alpha)][1:2] + +tsne.results <- matrix(rnorm(2*ncol(out1)), ncol=2) + +em1 <- embedVelocity(reducedDim(out1, 1), out1)[,1:2] +em2 <- embedVelocity(reducedDim(out2, 1), out2)[,1:2] + +out3 <- out1 +colData(out3) <- cbind(colData(out3), type = sample(letters, ncol(out3), replace = TRUE)) +rowData(out3) <- rowData(out3)[, -grep("gamma", colnames(rowData(out3)))] + +out4 <- out2 +rowData(out4) <- rowData(out4)[, -grep("fit_[us]0", colnames(rowData(out4)))] + +test_that("plotVelocity runs", { + + skip_if_not_installed("ggplot2") + skip_if_not_installed("patchwork") + + expect_error(plotVelocity("error", genes2)) + expect_error(plotVelocity(out1, "error")) + expect_error(plotVelocity(out1, genes2, use.dimred = "error")) + expect_error(plotVelocity(out1, genes2, use.dimred = FALSE)) + expect_error(plotVelocity(out1, genes2, 1, "error")) + expect_error(plotVelocity(out1, genes2, 1, which.plots = "error")) + expect_error(plotVelocity(out1, genes2, 1, genes.per.row = "error")) + expect_error(plotVelocity(out1, genes2, 1, color_by = c("error", "error"))) + expect_error(plotVelocity(out1, genes2, 1, max.abs.velo = "error")) + expect_error(plotVelocity(out1, genes2, 1, max.abs.velo = -1)) + + + tf <- tempfile(fileext = ".png") + png(tf) + expect_is(res0 <- plotVelocity(out1, genes1, which.plots = "phase", color_by = "#44444422"), "patchwork") + print(res0) + expect_is(res1 <- plotVelocity(out1, genes1), "patchwork") + print(res1) + expect_warning(res2 <- plotVelocity(out2, genes1)) + expect_is(res2, "patchwork") + print(res2) + expect_is(res3 <- plotVelocity(out4, genes2), "patchwork") + print(res3) + expect_is(res4 <- plotVelocity(out3, genes2, which.plots = "phase", color_by = "type"), "patchwork") + print(res4) + expect_is(res5 <- plotVelocity(out2, genes2), "patchwork") + print(res5) + dev.off() + expect_true(file.exists(tf)) + unlink(tf) +}) + +test_that("plotVelocityStream runs", { + + skip_if_not_installed("ggplot2") + skip_if_not_installed("metR") + + expect_error(plotVelocityStream("error", em2)) + expect_error(plotVelocityStream(out2, "error")) + expect_error(plotVelocityStream(out2, em2[1:10, ])) + expect_error(plotVelocityStream(out2, em2, use.dimred = "error")) + expect_error(plotVelocityStream(out2, em2, use.dimred = FALSE)) + expect_error(plotVelocityStream(out2, em2, color_by = "error")) + expect_error(plotVelocityStream(out2, em2, grid.resolution = "error")) + expect_error(plotVelocityStream(out2, em2, scale = "error")) + expect_error(plotVelocityStream(out2, em2, color.streamlines = "error")) + + tf <- tempfile(fileext = ".png") + png(tf) + expect_warning(print(plotVelocityStream(out2, em2, color_by = "#44444422"))) + print(plotVelocityStream(out2, em2)) + print(plotVelocityStream(out3, em2, color_by = "type")) + print(plotVelocityStream(out2, em2, color.streamlines = TRUE)) + dev.off() + expect_true(file.exists(tf)) + unlink(tf) +})