Multi-omics reveals heterogeneity and functional populations of oligodendrocyte progenitor cells induced by human neural stem cells.

Code ↔ Paper

5 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 5 matches

[1] § Materials and methods › Single-cell RNA sequencing (scRNA‑seq) of hOPCs ↔ R/visualization.R, lines 1035–1123 · score 0.76 · percent mito, genes detected, cells expressing, dimension reduction, mitochondrial, gene expression
[2] § Materials and methods › Data collection and analysis ↔ R/preprocessing.R, lines 2175–2251 · score 0.70 · Nuclear segmentation, DAPI staining, cytoplasmic, channel, intensity, exported
[3] § Materials and methods › Single-cell RNA sequencing (scRNA‑seq) of hOPCs ↔ R/doubletFinder.R, lines 1–66 · score 0.67 · Doublet Finder, principal components, Seurat, predicted, PCA, matrix
[4] § Materials and methods › GO and KEGG ↔ R/differential_expression.R, lines 455–538 · score 0.61 · DESeq2, fold change, v4, predict, threshold, gene
[5] § Materials and methods › GO and KEGG ↔ R/generics.R, lines 109–170 · score 0.57 · DESeq2, fold change, ClusterProfiler, PRE, gene

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

R · 6,161 lines · 195 KB · other · 1 match

#' @importFrom utils globalVariables
#' @importFrom ggplot2 fortify GeomViolin ggproto
#' @importFrom SeuratObject DefaultDimReduc
#'
NULL
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Generics
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#' @importFrom methods setGeneric
#'
setGeneric(
name = '.PrepImageData',
def = function(data, cells, ...) {
standardGeneric(f = '.PrepImageData')
}
)
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Heatmaps
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#' Dimensional reduction heatmap
#'
#' Draws a heatmap focusing on a principal component. Both cells and genes are sorted by their
#' principal component scores. Allows for nice visualization of sources of heterogeneity in the dataset.
#'
#' @inheritParams DoHeatmap
#' @param dims Dimensions to plot
#' @param nfeatures Number of genes to plot
#' @param cells A list of cells to plot. If numeric, just plots the top cells.
#' @param reduction Which dimensional reduction to use
#' @param balanced Plot an equal number of genes with both + and - scores.
#' @param projected Use the full projected dimensional reduction
#' @param ncol Number of columns to plot
#' @param fast If true, use \code{image} to generate plots; faster than using ggplot2,
#' but not customizable and excludes figure legend in output
#' @param assays A vector of assays to pull data from
#' @param combine Combine plots into a single \code{\link[patchwork]{patchwork}ed} ggplot object with
#' single shared figure legend when \code{fast=FALSE}. If \code{FALSE}, return a list of ggplot objects
#' @param legend.position When \code{combine=TRUE}, allows legend position to be adjusted
#' for \code{\link[patchwork]{patchwork}ed} output (default "right"). See \link[ggplot2]{theme}
#'
#' @return No return value by default. If using fast = FALSE, will return a
#' \code{\link[patchwork]{patchwork}ed} ggplot object if combine = TRUE, otherwise
#' returns a list of ggplot objects
#'
#' @importFrom patchwork wrap_plots plot_layout
#' @export
#' @concept visualization
#'
#' @seealso \code{\link[graphics]{image}} \code{\link[ggplot2]{geom_raster}}
#'
#' @examples
#' data("pbmc_small")
#' DimHeatmap(object = pbmc_small)
#'
DimHeatmap <- function(
object,
dims = 1,
nfeatures = 30,
cells = NULL,
reduction = 'pca',
disp.min = -2.5,
disp.max = NULL,
balanced = TRUE,
projected = FALSE,
ncol = NULL,
fast = TRUE,
raster = TRUE,
slot = 'scale.data',
assays = NULL,
combine = TRUE,
legend.position = "right"
) {
ncol <- ncol %||% ifelse(test = length(x = dims) > 2, yes = 3, no = length(x = dims))
plots <- vector(mode = 'list', length = length(x = dims))
assays <- assays %||% DefaultAssay(object = object)
disp.max <- disp.max %||% ifelse(
test = slot == 'scale.data',
yes = 2.5,
no = 6
)
if (!DefaultAssay(object = object[[reduction]]) %in% assays) {
warning("The original assay that the reduction was computed on is different than the assay specified")
}
cells <- cells %||% ncol(x = object)
if (is.numeric(x = cells)) {
cells <- lapply(
X = dims,
FUN = function(x) {
cells <- TopCells(
object = object[[reduction]],
dim = x,
ncells = cells,
balanced = balanced
)
if (balanced) {
cells$negative <- rev(x = cells$negative)
}
cells <- unlist(x = unname(obj = cells))
return(cells)
}
)
}
if (!is.list(x = cells)) {
cells <- lapply(X = 1:length(x = dims), FUN = function(x) {return(cells)})
}
features <- lapply(
X = dims,
FUN = TopFeatures,
object = object[[reduction]],
nfeatures = nfeatures,
balanced = balanced,
projected = projected
)
features.all <- unique(x = unlist(x = features))
if (length(x = assays) > 1) {
features.keyed <- lapply(
X = assays,
FUN = function(assay) {
features <- features.all[features.all %in% rownames(x = object[[assay]])]
if (length(x = features) > 0) {
return(paste0(Key(object = object[[assay]]), features))
}
}
)
features.keyed <- Filter(f = Negate(f = is.null), x = features.keyed)
features.keyed <- unlist(x = features.keyed)
} else {
features.keyed <- features.all
DefaultAssay(object = object) <- assays
}
data.all <- FetchData(
object = object,
vars = features.keyed,
cells = unique(x = unlist(x = cells)),
layer = slot
)
data.all <- MinMax(data = data.all, min = disp.min, max = disp.max)
data.limits <- c(min(data.all), max(data.all))
# if (check.plot && any(c(length(x = features.keyed), length(x = cells[[1]])) > 700)) {
# choice <- menu(c("Continue with plotting", "Quit"), title = "Plot(s) requested will likely take a while to plot.")
# if (choice != 1) {
# return(invisible(x = NULL))
# }
# }
if (fast) {
nrow <- floor(x = length(x = dims) / 3.01) + 1
orig.par <- par()$mfrow
par(mfrow = c(nrow, ncol))
}
for (i in 1:length(x = dims)) {
dim.features <- c(features[[i]][[2]], rev(x = features[[i]][[1]]))
dim.features <- rev(x = unlist(x = lapply(
X = dim.features,
FUN = function(feat) {
return(grep(pattern = paste0(feat, '$'), x = features.keyed, value = TRUE))
}
)))
dim.cells <- cells[[i]]
data.plot <- data.all[dim.cells, dim.features]
if (fast) {
SingleImageMap(
data = data.plot,
title = paste0(Key(object = object[[reduction]]), dims[i]),
order = dim.cells
)
} else {
plots[[i]] <- SingleRasterMap(
data = data.plot,
raster = raster,
limits = data.limits,
cell.order = dim.cells,
feature.order = dim.features
)
plots[[i]] <- plots[[i]] +
ggtitle(paste0(Key(object = object[[reduction]]), dims[i])) +
theme(plot.title = element_text(hjust = 0.5, face = "bold"))
}
}
if (fast) {
par(mfrow = orig.par)
return(invisible(x = NULL))
}
if (combine) {
plots <- wrap_plots(plots, ncol = ncol, guides = "collect") +
plot_layout(guides = "collect") &
theme(legend.position = legend.position)
}
return(plots)
}
#' Feature expression heatmap
#'
#' Draws a heatmap of single cell feature expression.
#'
#' @param object Seurat object
#' @param features A vector of features to plot, defaults to \code{VariableFeatures(object = object)}
#' @param cells A vector of cells to plot
#' @param disp.min Minimum display value (all values below are clipped)
#' @param disp.max Maximum display value (all values above are clipped); defaults to 2.5
#' if \code{slot} is 'scale.data', 6 otherwise
#' @param group.by A vector of variables to group cells by; pass 'ident' to group by cell identity classes
#' @param group.bar Add a color bar showing group status for cells
#' @param group.colors Colors to use for the color bar
#' @param slot Data slot to use, choose from 'raw.data', 'data', or 'scale.data'
#' @param assay Assay to pull from
# @param check.plot Check that plotting will finish in a reasonable amount of time
#' @param label Label the cell identies above the color bar
#' @param size Size of text above color bar
#' @param hjust Horizontal justification of text above color bar
#' @param vjust Vertical justification of text above color bar
#' @param angle Angle of text above color bar
#' @param raster If true, plot with geom_raster, else use geom_tile. geom_raster may look blurry on
#' some viewing applications such as Preview due to how the raster is interpolated. Set this to FALSE
#' if you are encountering that issue (note that plots may take longer to produce/render).
#' @param draw.lines Include white lines to separate the groups
#' @param lines.width Integer number to adjust the width of the separating white lines.
#' Corresponds to the number of "cells" between each group.
#' @param group.bar.height Scale the height of the color bar
#' @param combine Combine plots into a single \code{\link[patchwork]{patchwork}ed}
#' ggplot object. If \code{FALSE}, return a list of ggplot objects
#'
#' @return A \code{\link[patchwork]{patchwork}ed} ggplot object if
#' \code{combine = TRUE}; otherwise, a list of ggplot objects
#'
#' @importFrom stats median
#' @importFrom scales hue_pal
#' @importFrom ggplot2 annotation_raster coord_cartesian scale_color_manual
#' ggplot_build geom_text
#' @importFrom patchwork wrap_plots
#' @export
#' @concept visualization
#'
#' @examples
#' data("pbmc_small")
#' DoHeatmap(object = pbmc_small)
#'
DoHeatmap <- function(
object,
features = NULL,
cells = NULL,
group.by = 'ident',
group.bar = TRUE,
group.colors = NULL,
disp.min = -2.5,
disp.max = NULL,
slot = 'scale.data',
assay = NULL,
label = TRUE,
size = 5.5,
hjust = 0,
vjust = 0,
angle = 45,
raster = TRUE,
draw.lines = TRUE,
lines.width = NULL,
group.bar.height = 0.02,
combine = TRUE
) {
assay <- assay %||% DefaultAssay(object = object)
DefaultAssay(object = object) <- assay
cells <- cells %||% colnames(x = object[[assay]])
if (is.numeric(x = cells)) {
cells <- colnames(x = object)[cells]
}
features <- features %||% VariableFeatures(object = object)
features <- rev(x = unique(x = features))
disp.max <- disp.max %||% ifelse(
test = slot == 'scale.data',
yes = 2.5,
no = 6
)
# make sure features are present
possible.features <- Features(object,layer = slot)
if (any(!features %in% possible.features)) {
bad.features <- features[!features %in% possible.features]
features <- features[features %in% possible.features]
if(length(x = features) == 0) {
stop("No requested features found in the ", slot, " layer for the ", assay, " assay.")
}
warning("The following features were omitted as they were not found in the ", slot,
" layer for the ", assay, " assay: ", paste(bad.features, collapse = ", "))
}
data <- as.data.frame(x = as.matrix(x = t(x = GetAssayData(
object = object,
layer = slot)[features, cells, drop = FALSE])))
object <- suppressMessages(expr = StashIdent(object = object, save.name = 'ident'))
group.by <- group.by %||% 'ident'
groups.use <- object[[group.by]][cells, , drop = FALSE]
# group.use <- switch(
# EXPR = group.by,
# 'ident' = Idents(object = object),
# object[[group.by, drop = TRUE]]
# )
# group.use <- factor(x = group.use[cells])
plots <- vector(mode = 'list', length = ncol(x = groups.use))
for (i in 1:ncol(x = groups.use)) {
data.group <- data
group.use <- groups.use[, i, drop = TRUE]
if (!is.factor(x = group.use)) {
group.use <- factor(x = group.use)
}
names(x = group.use) <- cells
if (draw.lines) {
# create fake cells to serve as the white lines, fill with NAs
lines.width <- lines.width %||% ceiling(x = nrow(x = data.group) * 0.0025)
placeholder.cells <- sapply(
X = 1:(length(x = levels(x = group.use)) * lines.width),
FUN = function(x) {
return(RandomName(length = 20))
}
)
placeholder.groups <- rep(x = levels(x = group.use), times = lines.width)
group.levels <- levels(x = group.use)
names(x = placeholder.groups) <- placeholder.cells
group.use <- as.vector(x = group.use)
names(x = group.use) <- cells
group.use <- factor(x = c(group.use, placeholder.groups), levels = group.levels)
na.data.group <- matrix(
data = NA,
nrow = length(x = placeholder.cells),
ncol = ncol(x = data.group),
dimnames = list(placeholder.cells, colnames(x = data.group))
)
data.group <- rbind(data.group, na.data.group)
}
lgroup <- length(levels(group.use))
plot <- SingleRasterMap(
data = data.group,
raster = raster,
disp.min = disp.min,
disp.max = disp.max,
feature.order = features,
cell.order = names(x = sort(x = group.use)),
group.by = group.use
)
if (group.bar) {
# TODO: Change group.bar to annotation.bar
default.colors <- c(hue_pal()(length(x = levels(x = group.use))))
if (!is.null(x = names(x = group.colors))) {
cols <- unname(obj = group.colors[levels(x = group.use)])
} else {
cols <- group.colors[1:length(x = levels(x = group.use))] %||% default.colors
}
if (any(is.na(x = cols))) {
cols[is.na(x = cols)] <- default.colors[is.na(x = cols)]
cols <- Col2Hex(cols)
col.dups <- sort(x = unique(x = which(x = duplicated(x = substr(
x = cols,
start = 1,
stop = 7
)))))
through <- length(x = default.colors)
while (length(x = col.dups) > 0) {
pal.max <- length(x = col.dups) + through
cols.extra <- hue_pal()(pal.max)[(through + 1):pal.max]
cols[col.dups] <- cols.extra
col.dups <- sort(x = unique(x = which(x = duplicated(x = substr(
x = cols,
start = 1,
stop = 7
)))))
}
}
group.use2 <- sort(x = group.use)
if (draw.lines) {
na.group <- RandomName(length = 20)
levels(x = group.use2) <- c(levels(x = group.use2), na.group)
group.use2[placeholder.cells] <- na.group
cols <- c(cols, "#FFFFFF")
}
pbuild <- ggplot_build(plot = plot)
names(x = cols) <- levels(x = group.use2)
# scale the height of the bar
y.range <- diff(x = pbuild$layout$panel_params[[1]]$y.range)
y.pos <- max(pbuild$layout$panel_params[[1]]$y.range) + y.range * 0.015
y.max <- y.pos + group.bar.height * y.range
x.min <- min(pbuild$layout$panel_params[[1]]$x.range) + 0.1
x.max <- max(pbuild$layout$panel_params[[1]]$x.range) - 0.1
plot <- plot +
annotation_raster(
raster = t(x = cols[group.use2]),
xmin = x.min,
xmax = x.max,
ymin = y.pos,
ymax = y.max
) +
coord_cartesian(ylim = c(0, y.max), clip = 'off') +
scale_color_manual(
values = cols[-length(x = cols)],
name = "Identity",
na.translate = FALSE
)
if (label) {
x.max <- max(pbuild$layout$panel_params[[1]]$x.range)
# Attempt to pull xdivs from x.major in ggplot2 < 3.3.0; if NULL, pull from the >= 3.3.0 slot
x.divs <- pbuild$layout$panel_params[[1]]$x.major %||% attr(x = pbuild$layout$panel_params[[1]]$x$get_breaks(), which = "pos")
x <- data.frame(group = sort(x = group.use), x = x.divs)
label.x.pos <- tapply(X = x$x, INDEX = x$group, FUN = function(y) {
if (isTRUE(x = draw.lines)) {
mean(x = y[-length(x = y)])
} else {
mean(x = y)
}
})
label.x.pos <- data.frame(group = names(x = label.x.pos), label.x.pos)
plot <- plot + geom_text(
stat = "identity",
data = label.x.pos,
aes(label = .data[['group']], x = .data[['label.x.pos']]),
y = y.max + y.max * 0.03 * 0.5 + vjust,
angle = angle,
hjust = hjust,
size = size
)
plot <- suppressMessages(plot + coord_cartesian(
ylim = c(0, y.max + y.max * 0.002 * max(nchar(x = levels(x = group.use))) * size),
clip = 'off')
)
}
}
plot <- plot + theme(line = element_blank())
plots[[i]] <- plot
}
if (combine) {
plots <- wrap_plots(plots)
}
return(plots)
}
#' Hashtag oligo heatmap
#'
#' Draws a heatmap of hashtag oligo signals across singlets/doublets/negative cells. Allows for the visualization of HTO demultiplexing results.
#'
#' @param object Seurat object. Assumes that the hash tag oligo (HTO) data has been added and normalized, and demultiplexing has been run with HTODemux().
#' @param classification The naming for metadata column with classification result from HTODemux().
#' @param global.classification The slot for metadata column specifying a cell as singlet/doublet/negative.
#' @param assay Hashtag assay name.
#' @param ncells Number of cells to plot. Default is to choose 5000 cells by random subsampling, to avoid having to draw exceptionally large heatmaps.
#' @param singlet.names Namings for the singlets. Default is to use the same names as HTOs.
#' @param raster If true, plot with geom_raster, else use geom_tile. geom_raster may look blurry on
#' some viewing applications such as Preview due to how the raster is interpolated. Set this to FALSE
#' if you are encountering that issue (note that plots may take longer to produce/render).
#' @return Returns a ggplot2 plot object.
#'
#' @importFrom ggplot2 guides
#' @export
#' @concept visualization
#'
#' @seealso \code{\link{HTODemux}}
#'
#' @examples
#' \dontrun{
#' object <- HTODemux(object)
#' HTOHeatmap(object)
#' }
#'
HTOHeatmap <- function(
object,
assay = 'HTO',
classification = paste0(assay, '_classification'),
global.classification = paste0(assay, '_classification.global'),
ncells = 5000,
singlet.names = NULL,
raster = TRUE
) {
DefaultAssay(object = object) <- assay
Idents(object = object) <- object[[classification, drop = TRUE]]
if (ncells > ncol(x = object)) {
warning("ncells (", ncells, ") is larger than the number of cells present in the provided object (", ncol(x = object), "). Plotting heatmap for all cells.")
} else {
object <- subset(
x = object,
cells = sample(x = colnames(x = object), size = ncells)
)
}
classification <- object[[classification]]
singlets <- which(x = object[[global.classification]] == 'Singlet')
singlet.ids <- sort(x = unique(x = as.character(x = classification[singlets, ])))
doublets <- which(object[[global.classification]] == 'Doublet')
doublet.ids <- sort(x = unique(x = as.character(x = classification[doublets, ])))
heatmap.levels <- c(singlet.ids, doublet.ids, 'Negative')
object <- ScaleData(object = object, assay = assay, verbose = FALSE)
data <- FetchData(object = object, vars = singlet.ids)
Idents(object = object) <- factor(x = classification[, 1], levels = heatmap.levels)
plot <- SingleRasterMap(
data = data,
raster = raster,
feature.order = rev(x = singlet.ids),
cell.order = names(x = sort(x = Idents(object = object))),
group.by = Idents(object = object)
) + guides(color = "none")
return(plot)
}
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Expression by identity plots
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#' Single cell ridge plot
#'
#' Draws a ridge plot of single cell data (gene expression, metrics, PC
#' scores, etc.)
#'
#' @param object Seurat object
#' @param features Features to plot (gene expression, metrics, PC scores,
#' anything that can be retreived by FetchData)
#' @param cols Colors to use for plotting
#' @param idents Which classes to include in the plot (default is all)
#' @param sort Sort identity classes (on the x-axis) by the average
#' expression of the attribute being potted, can also pass 'increasing' or 'decreasing' to change sort direction
#' @param assay Name of assay to use, defaults to the active assay
#' @param group.by Group (color) cells in different ways (for example, orig.ident)
#' @param y.max Maximum y axis value
#' @param same.y.lims Set all the y-axis limits to the same values
#' @param log plot the feature axis on log scale
#' @param ncol Number of columns if multiple plots are displayed
#' @param slot Slot to pull expression data from (e.g. "counts" or "data")
#' @param layer Layer to pull expression data from (e.g. "counts" or "data")
#' @param stack Horizontally stack plots for each feature
#' @param combine Combine plots into a single \code{\link[patchwork]{patchwork}ed}
#' ggplot object. If \code{FALSE}, return a list of ggplot
#' @param fill.by Color violins/ridges based on either 'feature' or 'ident'
#'
#' @return A \code{\link[patchwork]{patchwork}ed} ggplot object if
#' \code{combine = TRUE}; otherwise, a list of ggplot objects
#'
#' @export
#' @concept visualization
#'
#' @examples
#' data("pbmc_small")
#' RidgePlot(object = pbmc_small, features = 'PC_1')
#'
RidgePlot <- function(
object,
features,
cols = NULL,
idents = NULL,
sort = FALSE,
assay = NULL,
group.by = NULL,
y.max = NULL,
same.y.lims = FALSE,
log = FALSE,
ncol = NULL,
slot = deprecated(),
layer = 'data',
stack = FALSE,
combine = TRUE,
fill.by = NULL
) {
if (is_present(arg = slot)) {
deprecate_soft(
when = '5.0.0',
what = 'RidgePlot(slot = )',
with = 'RidgePlot(layer = )'
)
layer <- slot %||% layer
}
fill.by <- fill.by %||% if (stack) 'feature' else 'ident'
return(ExIPlot(
object = object,
type = 'ridge',
features = features,
idents = idents,
ncol = ncol,
sort = sort,
assay = assay,
y.max = y.max,
same.y.lims = same.y.lims,
cols = cols,
group.by = group.by,
log = log,
layer = layer,
stack = stack,
combine = combine,
fill.by = fill.by
))
}
#' Single cell violin plot
#'
#' Draws a violin plot of single cell data (gene expression, metrics, PC
#' scores, etc.)
#'
#' @inheritParams RidgePlot
#' @param pt.size Point size for points
#' @param alpha Alpha value for points
#' @param split.by A factor in object metadata to split the plot by, pass 'ident'
#' to split by cell identity
#' @param split.plot plot each group of the split violin plots by multiple or
#' single violin shapes.
#' @param adjust Adjust parameter for geom_violin
#' @param flip flip plot orientation (identities on x-axis)
#' @param add.noise determine if adding a small noise for plotting
#' @param raster Convert points to raster format. Requires 'ggrastr' to be installed.
# default is \code{NULL} which automatically rasterizes if ggrastr is installed and
# number of points exceed 100,000.
#' @param raster.dpi the dpi for raster layer, default is 300.
#' See \code{\link[ggrastr]{rasterize}} for more info.
#'
#' @return A \code{\link[patchwork]{patchwork}ed} ggplot object if
#' \code{combine = TRUE}; otherwise, a list of ggplot objects
#'
#' @export
#' @concept visualization
#'
#' @seealso \code{\link{FetchData}}
#'
#' @examples
#' data("pbmc_small")
#' VlnPlot(object = pbmc_small, features = 'PC_1')
#' VlnPlot(object = pbmc_small, features = 'LYZ', split.by = 'groups')
#'
VlnPlot <- function(
object,
features,
cols = NULL,
pt.size = NULL,
alpha = 1,
idents = NULL,
sort = FALSE,
assay = NULL,
group.by = NULL,
split.by = NULL,
adjust = 1,
y.max = NULL,
same.y.lims = FALSE,
log = FALSE,
ncol = NULL,
slot = deprecated(),
layer = NULL,
split.plot = FALSE,
stack = FALSE,
combine = TRUE,
fill.by = 'feature',
flip = FALSE,
add.noise = TRUE,
raster = NULL,
raster.dpi = 300
) {
if (is_present(arg = slot)) {
deprecate_soft(
when = '5.0.0',
what = 'VlnPlot(slot = )',
with = 'VlnPlot(layer = )'
)
layer <- slot %||% layer
}
layer.set <- suppressWarnings(
Layers(
object = object,
search = layer %||% 'data'
)
)
if (is.null(layer) && length(layer.set) == 1 && layer.set == 'scale.data'){
warning('Default search for "data" layer yielded no results; utilizing "scale.data" layer instead.')
}
assay.name <- assay %||% DefaultAssay(object = object)
if (is.null(layer.set) & is.null(layer) ) {
warning('Default search for "data" layer in "', assay.name, '" assay yielded no results; utilizing "counts" layer instead.',
call. = FALSE, immediate. = TRUE)
layer.set <- Layers(
object = object,
search = 'counts'
)
}
if (is.null(layer.set)) {
stop('layer "', layer,'" is not found in assay: "', assay.name, '"')
} else {
layer <- layer.set
}
if (
!is.null(x = split.by) &
getOption(x = 'Seurat.warn.vlnplot.split', default = TRUE)
) {
message(
"The default behaviour of split.by has changed.\n",
"Separate violin plots are now plotted side-by-side.\n",
"To restore the old behaviour of a single split violin,\n",
"set split.plot = TRUE.
\nThis message will be shown once per session."
)
options(Seurat.warn.vlnplot.split = FALSE)
}
return(ExIPlot(
object = object,
type = ifelse(test = split.plot, yes = 'splitViolin', no = 'violin'),
features = features,
idents = idents,
ncol = ncol,
sort = sort,
assay = assay,
y.max = y.max,
same.y.lims = same.y.lims,
adjust = adjust,
pt.size = pt.size,
alpha = alpha,
cols = cols,
group.by = group.by,
split.by = split.by,
log = log,
layer = layer,
stack = stack,
combine = combine,
fill.by = fill.by,
flip = flip,
add.noise = add.noise,
raster = raster,
raster.dpi = raster.dpi
))
}
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Dimensional reduction plots
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#' Color dimensional reduction plot by tree split
#'
#' Returns a DimPlot colored based on whether the cells fall in clusters
#' to the left or to the right of a node split in the cluster tree.
#'
#' @param object Seurat object
#' @param node Node in cluster tree on which to base the split
#' @param left.color Color for the left side of the split
#' @param right.color Color for the right side of the split
#' @param other.color Color for all other cells
#' @inheritDotParams DimPlot -object
#'
#' @return Returns a DimPlot
#'
#' @export
#' @concept visualization
#'
#' @seealso \code{\link{DimPlot}}
#'
#' @examples
#' \dontrun{
#' if (requireNamespace("ape", quietly = TRUE)) {
#' data("pbmc_small")
#' pbmc_small <- BuildClusterTree(object = pbmc_small, verbose = FALSE)
#' PlotClusterTree(pbmc_small)
#' ColorDimSplit(pbmc_small, node = 5)
#' }
#' }
#'
ColorDimSplit <- function(
object,
node,
left.color = 'red',
right.color = 'blue',
other.color = 'grey50',
...
) {
CheckDots(..., fxns = 'DimPlot')
tree <- Tool(object = object, slot = "BuildClusterTree")
split <- tree$edge[which(x = tree$edge[, 1] == node), ][, 2]
all.children <- sort(x = tree$edge[, 2][!tree$edge[, 2] %in% tree$edge[, 1]])
left.group <- DFT(tree = tree, node = split[1], only.children = TRUE)
right.group <- DFT(tree = tree, node = split[2], only.children = TRUE)
if (any(is.na(x = left.group))) {
left.group <- split[1]
}
if (any(is.na(x = right.group))) {
right.group <- split[2]
}
left.group <- MapVals(v = left.group, from = all.children, to = tree$tip.label)
right.group <- MapVals(v = right.group, from = all.children, to = tree$tip.label)
remaining.group <- setdiff(x = tree$tip.label, y = c(left.group, right.group))
left.cells <- WhichCells(object = object, ident = left.group)
right.cells <- WhichCells(object = object, ident = right.group)
remaining.cells <- WhichCells(object = object, ident = remaining.group)
object <- SetIdent(
object = object,
cells = left.cells,
value = "Left Split"
)
object <- SetIdent(
object = object,
cells = right.cells,
value = "Right Split"
)
object <- SetIdent(
object = object,
cells = remaining.cells,
value = "Not in Split"
)
levels(x = object) <- c("Left Split", "Right Split", "Not in Split")
colors.use = c(left.color, right.color, other.color)
return(DimPlot(object = object, cols = colors.use, ...))
}
#' Dimensional reduction plot
#'
#' Graphs the output of a dimensional reduction technique on a 2D scatter plot where each point is a
#' cell and it's positioned based on the cell embeddings determined by the reduction technique. By
#' default, cells are colored by their identity class (can be changed with the group.by parameter).
#'
#' @param object Seurat object
#' @param dims Dimensions to plot, must be a two-length numeric vector specifying x- and y-dimensions
#' @param cells Vector of cells to plot (default is all cells)
#' @param cols Vector of colors, each color corresponds to an identity class. This may also be a single character
#' or numeric value corresponding to a palette as specified by \code{\link[RColorBrewer]{brewer.pal.info}}.
#' By default, ggplot2 assigns colors. We also include a number of palettes from the pals package.
#' See \code{\link{DiscretePalette}} for details.
#' @param pt.size Adjust point size for plotting
#' @param reduction Which dimensionality reduction to use. If not specified, first searches for umap, then tsne, then pca
#' @param group.by Name of one or more metadata columns to group (color) cells by
#' (for example, orig.ident); pass 'ident' to group by identity class
#' @param split.by A factor in object metadata to split the plot by, pass 'ident'
#' to split by cell identity
#' @param shape.by If NULL, all points are circles (default). You can specify any
#' cell attribute (that can be pulled with FetchData) allowing for both
#' different colors and different shapes on cells. Only applicable if \code{raster = FALSE}.
#' @param order Specify the order of plotting for the idents. This can be
#' useful for crowded plots if points of interest are being buried. Provide
#' either a full list of valid idents or a subset to be plotted last (on top)
#' @param shuffle Whether to randomly shuffle the order of points. This can be
#' useful for crowded plots if points of interest are being buried. (default is FALSE)
#' @param seed Sets the seed if randomly shuffling the order of points.
#' @param label Whether to label the clusters
#' @param label.size Sets size of labels
#' @param label.color Sets the color of the label text
#' @param label.box Whether to put a box around the label text (geom_text vs
#' geom_label)
#' @param alpha Alpha value for plotting (default is 1)
#' @param repel Repel labels
#' @param stroke.size Adjust stroke (outline) size of points
#' @param cells.highlight A list of character or numeric vectors of cells to
#' highlight. If only one group of cells desired, can simply
#' pass a vector instead of a list. If set, colors selected cells to the color(s)
#' in \code{cols.highlight} and other cells black (white if dark.theme = TRUE);
#' will also resize to the size(s) passed to \code{sizes.highlight}
#' @param cols.highlight A vector of colors to highlight the cells as; will
#' repeat to the length groups in cells.highlight
#' @param sizes.highlight Size of highlighted cells; will repeat to the length
#' groups in cells.highlight. If \code{sizes.highlight = TRUE} size of all
#' points will be this value.
#' @param na.value Color value for NA points when using custom scale
#' @param ncol Number of columns for display when combining plots
#' @param combine Combine plots into a single \code{\link[patchwork]{patchwork}ed}
#' ggplot object. If \code{FALSE}, return a list of ggplot objects
#' @param raster Convert points to raster format, default is \code{NULL} which
#' automatically rasterizes if plotting more than 100,000 cells
#' @param raster.dpi Pixel resolution for rasterized plots, passed to geom_scattermore(). Default is c(512, 512).
#' @param label.size.cutoff Clusters with fewer cells than the cutoff are not labeled (replaced with ' ' label)
#'
#' @return A \code{\link[patchwork]{patchwork}ed} ggplot object if
#' \code{combine = TRUE}; otherwise, a list of ggplot objects
#'
#' @importFrom rlang !!
#' @importFrom ggplot2 facet_wrap vars sym labs
#' @importFrom patchwork wrap_plots
#'
#' @export
#' @concept visualization
#'
#' @note For the old \code{do.hover} and \code{do.identify} functionality, please see
#' \code{HoverLocator} and \code{CellSelector}, respectively.
#'
#' @aliases TSNEPlot PCAPlot ICAPlot
#' @seealso \code{\link{FeaturePlot}} \code{\link{HoverLocator}}
#' \code{\link{CellSelector}} \code{\link{FetchData}}
#'
#' @examples
#' data("pbmc_small")
#' DimPlot(object = pbmc_small)
#' DimPlot(object = pbmc_small, split.by = 'letter.idents')
#'
DimPlot <- function(
object,
dims = c(1, 2),
cells = NULL,
cols = NULL,
pt.size = NULL,
reduction = NULL,
group.by = NULL,
split.by = NULL,
shape.by = NULL,
order = NULL,
shuffle = FALSE,
seed = 1,
label = FALSE,
label.size = 4,
label.color = 'black',
label.box = FALSE,
repel = FALSE,
alpha = 1,
stroke.size = NULL,
cells.highlight = NULL,
cols.highlight = '#DE2D26',
sizes.highlight = 1,
na.value = 'grey50',
ncol = NULL,
combine = TRUE,
raster = NULL,
raster.dpi = c(512, 512),
label.size.cutoff = 0
) {
if (!is_integerish(x = dims, n = 2L, finite = TRUE) || !all(dims > 0L)) {
abort(message = "'dims' must be a two-length integer vector")
}
reduction <- reduction %||% DefaultDimReduc(object = object)
# cells <- cells %||% colnames(x = object)
##### Cells for all cells in the assay.
#### Cells function should not only get default layer
cells <- cells %||% Cells(
x = object,
assay = DefaultAssay(object = object[[reduction]])
)
# data <- Embeddings(object = object[[reduction]])[cells, dims]
# data <- as.data.frame(x = data)
dims <- paste0(Key(object = object[[reduction]]), dims)
orig.groups <- group.by
group.by <- group.by %||% 'ident'
if (label & (label.size.cutoff > 0)) {
labels <- FetchData(object, group.by)
for(i in seq_along(group.by)) {
grouping_var <- group.by[i]
label_table <- table(labels[,grouping_var])
invalid_labels <- names(which(label_table < label.size.cutoff))
labels[, i] <- as.character(labels[,i])
labels[which(labels[,grouping_var] %in% invalid_labels),grouping_var] <- ' '
colnames(labels)[i] <- paste0(colnames(labels)[i], "_filtered")
}
object <- AddMetaData(object,labels)
group.by <- colnames(labels)
}
# check for overlap between colnames of dim reduc embeddings and metadata
metadata_cols <- names(object[[]])
colname_overlap <- intersect(dims, metadata_cols)
if (length(colname_overlap) > 0) {
warning("Found metadata columns with the same names as requested reduction columns: ",
paste(colname_overlap, collapse = ", "),
". Consider renaming these metadata column(s) to avoid conflicts with dimensionality reduction embeddings.",
call. = FALSE)
}
data <- FetchData(
object = object,
vars = c(dims, group.by),
cells = cells,
clean = 'project'
)
# cells <- rownames(x = object)
# object[['ident']] <- Idents(object = object)
# orig.groups <- group.by
# group.by <- group.by %||% 'ident'
# data <- cbind(data, object[[group.by]][cells, , drop = FALSE])
group.by <- colnames(x = data)[3:ncol(x = data)]
for (group in group.by) {
if (!is.factor(x = data[, group])) {
data[, group] <- factor(x = data[, group])
}
}
if (!is.null(x = shape.by)) {
data[, shape.by] <- object[[shape.by, drop = TRUE]]
}
if (!is.null(x = split.by)) {
split <- FetchData(object = object, vars = split.by, clean=TRUE)[split.by]
data <- data[rownames(split),]
data[, split.by] <- split
}
if (isTRUE(x = shuffle)) {
set.seed(seed = seed)
data <- data[sample(x = 1:nrow(x = data)), ]
}
plots <- lapply(
X = group.by,
FUN = function(x) {
plot <- SingleDimPlot(
data = data[, c(dims, x, split.by, shape.by)],
dims = dims,
col.by = x,
cols = cols,
pt.size = pt.size,
shape.by = shape.by,
order = order,
alpha = alpha,
stroke.size = stroke.size,
label = FALSE,
cells.highlight = cells.highlight,
cols.highlight = cols.highlight,
sizes.highlight = sizes.highlight,
na.value = na.value,
raster = raster,
raster.dpi = raster.dpi
)
if (label) {
plot <- LabelClusters(
plot = plot,
id = x,
repel = repel,
size = label.size,
split.by = split.by,
box = label.box,
color = label.color
)
}
if (!is.null(x = split.by)) {
plot <- plot + FacetTheme() +
facet_wrap(
facets = vars(!!sym(x = split.by)),
ncol = if (length(x = group.by) > 1 || is.null(x = ncol)) {
length(x = unique(x = data[, split.by]))
} else {
ncol
}
)
}
plot <- if (is.null(x = orig.groups)) {
plot + labs(title = NULL)
} else {
plot + CenterTitle()
}
}
)
if (!is.null(x = split.by)) {
ncol <- 1
}
if (combine) {
plots <- wrap_plots(plots, ncol = orig.groups %iff% ncol)
}
return(plots)
}
#' Visualize 'features' on a dimensional reduction plot
#'
#' Colors single cells on a dimensional reduction plot according to a 'feature'
#' (i.e. gene expression, PC scores, number of genes detected, etc.)
#'
#' @inheritParams DimPlot
#' @param order Boolean determining whether to plot cells in order of expression. Can be useful if
#' cells expressing given feature are getting buried.
#' @param features Vector of features to plot. Features can come from:
#' \itemize{
#' \item An \code{Assay} feature (e.g. a gene name - "MS4A1")
#' \item A column name from meta.data (e.g. mitochondrial percentage -
#' "percent.mito")
#' \item A column name from a \code{DimReduc} object corresponding to the
#' cell embedding values (e.g. the PC 1 scores - "PC_1")
#' }
#' @param cols The two colors to form the gradient over. Provide as string vector with
#' the first color corresponding to low values, the second to high. Also accepts a Brewer
#' color scale or vector of colors. Note: this will bin the data into number of colors provided.
#' When blend is \code{TRUE}, takes anywhere from 1-3 colors:
#' \describe{
#' \item{1 color:}{Treated as color for double-negatives, will use default colors 2 and 3 for per-feature expression}
#' \item{2 colors:}{Treated as colors for per-feature expression, will use default color 1 for double-negatives}
#' \item{3+ colors:}{First color used for double-negatives, colors 2 and 3 used for per-feature expression, all others ignored}
#' }
#' @param min.cutoff,max.cutoff Vector of minimum and maximum cutoff values for each feature,
#' may specify quantile in the form of 'q##' where '##' is the quantile (eg, 'q1', 'q10')
#' @param stroke.size Adjust stroke (outline) size of points
#' @param split.by A factor in object metadata to split the plot by, pass 'ident'
#' to split by cell identity
#' @param keep.scale How to handle the color scale across multiple plots. Options are:
#' \itemize{
#' \item \dQuote{feature} (default; by row/feature scaling): The plots for
#' each individual feature are scaled to the maximum expression of the
#' feature across the conditions provided to \code{split.by}
#' \item \dQuote{all} (universal scaling): The plots for all features and
#' conditions are scaled to the maximum expression value for the feature
#' with the highest overall expression
#' \item \code{NULL} (no scaling): Each individual plot is scaled to the
#' maximum expression value of the feature in the condition provided to
#' \code{split.by}. Be aware setting \code{NULL} will result in color
#' scales that are not comparable between plots
#' }
#' @param slot Which slot to pull expression data from?
#' @param assay Primary assay to pull feature data from
#' @param blend Scale and blend expression values to visualize coexpression of two features
#' @param blend.threshold The color cutoff from weak signal to strong signal; ranges from 0 to 1.
#' @param ncol Number of columns to combine multiple feature plots to, ignored if \code{split.by} is not \code{NULL}
#' @param coord.fixed Plot cartesian coordinates with fixed aspect ratio
#' @param by.col If splitting by a factor, plot the splits per column with the features as rows; ignored if \code{blend = TRUE}
#' @param sort.cell Redundant with \code{order}. This argument is being
#' deprecated. Please use \code{order} instead.
#' @param interactive Launch an interactive \code{\link[Seurat:IFeaturePlot]{FeaturePlot}}
#' @param combine Combine plots into a single \code{\link[patchwork]{patchwork}ed}
#' ggplot object. If \code{FALSE}, return a list of ggplot objects
#'
#' @return A \code{\link[patchwork]{patchwork}ed} ggplot object if
#' \code{combine = TRUE}; otherwise, a list of ggplot objects
#'
#' @importFrom grDevices rgb
#' @importFrom patchwork wrap_plots
#' @importFrom cowplot theme_cowplot
#' @importFrom RColorBrewer brewer.pal.info
#' @importFrom ggplot2 labs scale_x_continuous scale_y_continuous theme element_rect
#' dup_axis guides element_blank element_text margin scale_color_brewer scale_color_gradientn
#' scale_color_manual coord_fixed ggtitle
#'
#' @export
#' @concept visualization
#'
#' @note For the old \code{do.hover} and \code{do.identify} functionality, please see
#' \code{HoverLocator} and \code{CellSelector}, respectively.
#'
#' @aliases FeatureHeatmap
#' @seealso \code{\link{DimPlot}} \code{\link{HoverLocator}}
#' \code{\link{CellSelector}}
#'
#' @examples
#' data("pbmc_small")
#' FeaturePlot(object = pbmc_small, features = 'PC_1')
#'
FeaturePlot <- function(
object,
features,
dims = c(1, 2),
cells = NULL,
cols = if (blend) {
c('lightgrey', '#ff0000', '#00ff00')
} else {
c('lightgrey', 'blue')
},
pt.size = NULL,
alpha = 1,
stroke.size = NULL,
order = FALSE,
min.cutoff = NA,
max.cutoff = NA,
reduction = NULL,
split.by = NULL,
keep.scale = "feature",
shape.by = NULL,
slot = 'data',
assay = NULL,
blend = FALSE,
blend.threshold = 0.5,
label = FALSE,
label.size = 4,
label.color = "black",
repel = FALSE,
ncol = NULL,
coord.fixed = FALSE,
by.col = TRUE,
sort.cell = deprecated(),
interactive = FALSE,
combine = TRUE,
raster = NULL,
raster.dpi = c(512, 512)
) {
# TODO: deprecate fully on 3.2.0
if (is_present(arg = sort.cell)) {
deprecate_stop(
when = '4.9.0',
what = 'FeaturePlot(sort.cell = )',
with = 'FeaturePlot(order = )'
)
}
if (isTRUE(x = interactive)) {
return(IFeaturePlot(
object = object,
feature = features[1],
dims = dims,
reduction = reduction,
slot = slot
))
}
# Check keep.scale param for valid entries
if (!is.null(x = keep.scale)) {
keep.scale <- arg_match0(arg = keep.scale, values = c('feature', 'all'))
}
# Set a theme to remove right-hand Y axis lines
# Also sets right-hand Y axis text label formatting
no.right <- theme(
axis.line.y.right = element_blank(),
axis.ticks.y.right = element_blank(),
axis.text.y.right = element_blank(),
axis.title.y.right = element_text(
face = "bold",
size = 14,
margin = margin(r = 7)
)
)
# Get the DimReduc to use
reduction <- reduction %||% DefaultDimReduc(object = object)
if (!is_integerish(x = dims, n = 2L, finite = TRUE) && !all(dims > 0L)) {
abort(message = "'dims' must be a two-length integer vector")
}
# Figure out blending stuff
if (isTRUE(x = blend) && length(x = features) != 2) {
abort(message = "Blending feature plots only works with two features")
}
# Set color scheme for blended FeaturePlots
if (isTRUE(x = blend)) {
default.colors <- eval(expr = formals(fun = FeaturePlot)$cols)
cols <- switch(
EXPR = as.character(x = length(x = cols)),
'0' = {
warn(message = "No colors provided, using default colors")
default.colors
},
'1' = {
warn(message = paste(
"Only one color provided, assuming",
sQuote(x = cols),
"is double-negative and augmenting with default colors"
))
c(cols, default.colors[2:3])
},
'2' = {
warn(message = paste(
"Only two colors provided, assuming specified are for features and agumenting with",
sQuote(default.colors[1]),
"for double-negatives",
))
c(default.colors[1], cols)
},
'3' = cols,
{
warn(message = "More than three colors provided, using only first three")
cols[1:3]
}
)
}
if (isTRUE(x = blend) && length(x = cols) != 3) {
abort("Blending feature plots only works with three colors; first one for negative cells")
}
# Name the reductions
dims <- paste0(Key(object = object[[reduction]]), dims)
cells <- cells %||% Cells(x = object[[reduction]])
# Get plotting data
data <- FetchData(
object = object,
vars = c(dims, 'ident', features),
cells = cells,
layer = slot,
assay = assay
)
# Check presence of features/dimensions
if (ncol(x = data) < 4) {
abort(message = paste(
"None of the requested features were found:",
paste(features, collapse = ', '),
"in slot ",
slot
))
} else if (!all(dims %in% colnames(x = data))) {
abort(message = "The dimensions requested were not found")
}
features <- setdiff(x = names(x = data), y = c(dims, 'ident'))
# Determine cutoffs
min.cutoff <- mapply(
FUN = function(cutoff, feature) {
return(ifelse(
test = is.na(x = cutoff),
yes = min(data[, feature]),
no = cutoff
))
},
cutoff = min.cutoff,
feature = features
)
max.cutoff <- mapply(
FUN = function(cutoff, feature) {
return(ifelse(
test = is.na(x = cutoff),
yes = max(data[, feature]),
no = cutoff
))
},
cutoff = max.cutoff,
feature = features
)
check.lengths <- unique(x = vapply(
X = list(features, min.cutoff, max.cutoff),
FUN = length,
FUN.VALUE = numeric(length = 1)
))
if (length(x = check.lengths) != 1) {
abort(
message = "There must be the same number of minimum and maximum cuttoffs as there are features"
)
}
names(x = min.cutoff) <- names(x = max.cutoff) <- features
brewer.gran <- ifelse(
test = length(x = cols) == 1,
yes = brewer.pal.info[cols, ]$maxcolors,
no = length(x = cols)
)
# Apply cutoffs
for (i in seq_along(along.with = features)) {
f <- features[i]
data.feature <- data[[f]]
min.use <- SetQuantile(cutoff = min.cutoff[f], data = data.feature)
max.use <- SetQuantile(cutoff = max.cutoff[f], data = data.feature)
data.feature[data.feature < min.use] <- min.use
data.feature[data.feature > max.use] <- max.use
if (brewer.gran != 2) {
data.feature <- if (all(data.feature == 0)) {
rep_len(x = 0, length.out = length(x = data.feature))
} else {
as.numeric(x = as.factor(x = cut(
x = as.numeric(x = data.feature),
breaks = 2
)))
}
}
data[[f]] <- data.feature
}
# Figure out splits (FeatureHeatmap)
data$split <- if (is.null(x = split.by)) {
RandomName()
} else {
switch(
EXPR = split.by,
ident = Idents(object = object)[cells, drop = TRUE],
object[[split.by, drop = TRUE]][cells, drop = TRUE]
)
}
if (!is.factor(x = data$split)) {
data$split <- factor(x = data$split)
}
# Set shaping variable
if (!is.null(x = shape.by)) {
data[, shape.by] <- object[[shape.by, drop = TRUE]]
}
# Make list of plots
plots <- vector(
mode = "list",
length = ifelse(
test = blend,
yes = 4,
no = length(x = features) * length(x = levels(x = data$split))
)
)
# Apply common limits
xlims <- c(floor(x = min(data[, dims[1]])), ceiling(x = max(data[, dims[1]])))
ylims <- c(floor(min(data[, dims[2]])), ceiling(x = max(data[, dims[2]])))
# Set blended colors
if (blend) {
ncol <- 4
color.matrix <- BlendMatrix(
two.colors = cols[2:3],
col.threshold = blend.threshold,
negative.color = cols[1]
)
cols <- cols[2:3]
colors <- list(
color.matrix[, 1],
color.matrix[1, ],
as.vector(x = color.matrix)
)
}
# Make the plots
for (i in 1:length(x = levels(x = data$split))) {
# Figure out which split we're working with
ident <- levels(x = data$split)[i]
data.plot <- data[as.character(x = data$split) == ident, , drop = FALSE]
# Blend expression values
if (isTRUE(x = blend)) {
features <- features[1:2]
no.expression <- features[colMeans(x = data.plot[, features]) == 0]
if (length(x = no.expression) != 0) {
abort(message = paste(
"The following features have no value:",
paste(no.expression, collapse = ', ')
))
}
data.plot <- cbind(data.plot[, c(dims, 'ident')], BlendExpression(data = data.plot[, features[1:2]]))
features <- colnames(x = data.plot)[4:ncol(x = data.plot)]
}
# Make per-feature plots
for (j in 1:length(x = features)) {
feature <- features[j]
# Get blended colors
if (isTRUE(x = blend)) {
cols.use <- as.numeric(x = as.character(x = data.plot[, feature])) + 1
cols.use <- colors[[j]][sort(x = unique(x = cols.use))]
} else {
cols.use <- NULL
}
data.single <- data.plot[, c(dims, 'ident', feature, shape.by)]
# Make the plot
plot <- SingleDimPlot(
data = data.single,
dims = dims,
col.by = feature,
order = order,
pt.size = pt.size,
alpha = alpha,
stroke.size = stroke.size,
cols = cols.use,
shape.by = shape.by,
label = FALSE,
raster = raster,
raster.dpi = raster.dpi
) +
scale_x_continuous(limits = xlims) +
scale_y_continuous(limits = ylims) +
theme_cowplot() +
CenterTitle()
# theme(plot.title = element_text(hjust = 0.5))
# Add labels
if (isTRUE(x = label)) {
plot <- LabelClusters(
plot = plot,
id = 'ident',
repel = repel,
size = label.size,
color = label.color
)
}
# Make FeatureHeatmaps look nice(ish)
if (length(x = levels(x = data$split)) > 1) {
plot <- plot + theme(panel.border = element_rect(fill = NA, colour = 'black'))
# Add title
plot <- plot + if (i == 1) {
labs(title = feature)
} else {
labs(title = NULL)
}
# Add second axis
if (j == length(x = features) && !blend) {
suppressMessages(
expr = plot <- plot +
scale_y_continuous(
sec.axis = dup_axis(name = ident),
limits = ylims
) +
no.right
)
}
# Remove left Y axis
if (j != 1) {
plot <- plot + theme(
axis.line.y = element_blank(),
axis.ticks.y = element_blank(),
axis.text.y = element_blank(),
axis.title.y.left = element_blank()
)
}
# Remove bottom X axis
if (i != length(x = levels(x = data$split))) {
plot <- plot + theme(
axis.line.x = element_blank(),
axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
axis.title.x = element_blank()
)
}
} else {
plot <- plot + labs(title = feature)
}
# Add colors scale for normal FeaturePlots
if (!blend) {
plot <- plot + guides(color = NULL)
cols.grad <- cols
if (length(x = cols) == 1) {
plot <- plot + scale_color_brewer(palette = cols)
} else if (length(x = cols) > 1) {
unique.feature.exp <- unique(data.plot[, feature])
if (length(unique.feature.exp) == 1) {
warn(message = paste0(
"All cells have the same value (",
unique.feature.exp,
") of ",
dQuote(x = feature)
))
if (unique.feature.exp == 0) {
cols.grad <- cols[1]
} else{
cols.grad <- cols
}
}
plot <- suppressMessages(
expr = plot + scale_color_gradientn(
colors = cols.grad,
guide = "colorbar"
)
)
}
}
if (!(is.null(x = keep.scale)) && keep.scale == "feature" && !blend) {
max.feature.value <- max(data[, feature])
min.feature.value <- min(data[, feature])
plot <- suppressMessages(plot & scale_color_gradientn(colors = cols, limits = c(min.feature.value, max.feature.value)))
}
# Add coord_fixed
if (coord.fixed) {
plot <- plot + coord_fixed()
}
# I'm not sure why, but sometimes the damn thing fails without this
# Thanks ggplot2
plot <- plot
# Place the plot
plots[[(length(x = features) * (i - 1)) + j]] <- plot
}
}
# Add blended color key
if (isTRUE(x = blend)) {
blend.legend <- BlendMap(color.matrix = color.matrix)
for (ii in 1:length(x = levels(x = data$split))) {
suppressMessages(expr = plots <- append(
x = plots,
values = list(
blend.legend +
scale_y_continuous(
sec.axis = dup_axis(name = ifelse(
test = length(x = levels(x = data$split)) > 1,
yes = levels(x = data$split)[ii],
no = ''
)),
expand = c(0, 0)
) +
labs(
x = features[1],
y = features[2],
title = if (ii == 1) {
paste('Color threshold:', blend.threshold)
} else {
NULL
}
) +
no.right
),
after = 4 * ii - 1
))
}
}
# Remove NULL plots
plots <- Filter(f = Negate(f = is.null), x = plots)
# Combine the plots
if (is.null(x = ncol)) {
ncol <- 2
if (length(x = features) == 1) {
ncol <- 1
}
if (length(x = features) > 6) {
ncol <- 3
}
if (length(x = features) > 9) {
ncol <- 4
}
}
ncol <- ifelse(
test = is.null(x = split.by) || isTRUE(x = blend),
yes = ncol,
no = length(x = features)
)
legend <- if (isTRUE(x = blend)) {
'none'
} else {
split.by %iff% 'none'
}
# Transpose the FeatureHeatmap matrix (not applicable for blended FeaturePlots)
if (isTRUE(x = combine)) {
if (by.col && !is.null(x = split.by) && !blend) {
plots <- lapply(
X = plots,
FUN = function(x) {
return(suppressMessages(
expr = x +
theme_cowplot() +
ggtitle("") +
scale_y_continuous(sec.axis = dup_axis(name = ""), limits = ylims) +
no.right
))
}
)
nsplits <- length(x = levels(x = data$split))
idx <- 1
for (i in (length(x = features) * (nsplits - 1) + 1):(length(x = features) * nsplits)) {
plots[[i]] <- suppressMessages(
expr = plots[[i]] +
scale_y_continuous(
sec.axis = dup_axis(name = features[[idx]]),
limits = ylims
) +
no.right
)
idx <- idx + 1
}
idx <- 1
for (i in which(x = 1:length(x = plots) %% length(x = features) == 1)) {
plots[[i]] <- plots[[i]] +
ggtitle(levels(x = data$split)[[idx]]) +
theme(plot.title = element_text(hjust = 0.5))
idx <- idx + 1
}
idx <- 1
if (length(x = features) == 1) {
for (i in 1:length(x = plots)) {
plots[[i]] <- plots[[i]] +
ggtitle(levels(x = data$split)[[idx]]) +
theme(plot.title = element_text(hjust = 0.5))
idx <- idx + 1
}
ncol <- 1
nrow <- nsplits
} else {
nrow <- split.by %iff% length(x = levels(x = data$split))
}
plots <- plots[c(do.call(
what = rbind,
args = split(
x = 1:length(x = plots),
f = ceiling(x = seq_along(along.with = 1:length(x = plots)) / length(x = features))
)
))]
# Set ncol to number of splits (nrow) and nrow to number of features (ncol)
plots <- wrap_plots(plots, ncol = nrow, nrow = ncol)
if (!is.null(x = legend) && legend == 'none') {
plots <- plots & NoLegend()
}
} else {
plots <- wrap_plots(plots, ncol = ncol, nrow = split.by %iff% length(x = levels(x = data$split)))
}
if (!is.null(x = legend) && legend == 'none') {
plots <- plots & NoLegend()
}
if (!(is.null(x = keep.scale)) && keep.scale == "all" && !blend) {
max.feature.value <- max(data[, features])
min.feature.value <- min(data[, features])
plots <- suppressMessages(plots & scale_color_gradientn(colors = cols, limits = c(min.feature.value, max.feature.value)))
}
}
return(plots)
}
#' Visualize features in dimensional reduction space interactively
#'
#' @inheritParams FeaturePlot
#' @param feature Feature to plot
#'
#' @return Returns the final plot as a ggplot object
#'
#' @importFrom cowplot theme_cowplot
#' @importFrom ggplot2 theme element_text guides scale_color_gradientn
#' @importFrom miniUI miniPage miniButtonBlock miniTitleBarButton miniContentPanel
#' @importFrom shiny fillRow sidebarPanel selectInput plotOutput reactiveValues
#' observeEvent stopApp observe updateSelectInput renderPlot runGadget
#'
#' @export
#' @concept visualization
#'
IFeaturePlot <- function(object, feature, dims = c(1, 2), reduction = NULL, slot = 'data') {
# Set initial data values
feature.label <- 'Feature to visualize'
assay.keys <- Key(object = object)[Assays(object = object)]
keyed <- sapply(X = assay.keys, FUN = grepl, x = feature)
assay <- if (any(keyed)) {
names(x = which(x = keyed))[1]
} else {
DefaultAssay(object = object)
}
features <- sort(x = rownames(x = GetAssayData(
object = object,
layer = slot,
assay = assay
)))
assays.use <- vapply(
X = Assays(object = object),
FUN = function(x) {
return(!IsMatrixEmpty(x = GetAssayData(
object = object,
layer = slot,
assay = x
)))
},
FUN.VALUE = logical(length = 1L)
)
assays.use <- sort(x = Assays(object = object)[assays.use])
reduction <- reduction %||% DefaultDimReduc(object = object)
dims.reduc <- gsub(
pattern = Key(object = object[[reduction]]),
replacement = '',
x = colnames(x = object[[reduction]])
)
# Set up the gadget UI
ui <- miniPage(
miniButtonBlock(miniTitleBarButton(
inputId = 'done',
label = 'Done',
primary = TRUE
)),
miniContentPanel(
fillRow(
sidebarPanel(
selectInput(
inputId = 'assay',
label = 'Assay',
choices = assays.use,
selected = assay,
selectize = FALSE,
width = '100%'
),
selectInput(
inputId = 'feature',
label = feature.label,
choices = features,
selected = feature,
selectize = FALSE,
width = '100%'
),
selectInput(
inputId = 'reduction',
label = 'Dimensional reduction',
choices = Reductions(object = object),
selected = reduction,
selectize = FALSE,
width = '100%'
),
selectInput(
inputId = 'xdim',
label = 'X dimension',
choices = dims.reduc,
selected = as.character(x = dims[1]),
selectize = FALSE,
width = '100%'
),
selectInput(
inputId = 'ydim',
label = 'Y dimension',
choices = dims.reduc,
selected = as.character(x = dims[2]),
selectize = FALSE,
width = '100%'
),
selectInput(
inputId = 'palette',
label = 'Color scheme',
choices = names(x = FeaturePalettes),
selected = 'Seurat',
selectize = FALSE,
width = '100%'
),
width = '100%'
),
plotOutput(outputId = 'plot', height = '100%'),
flex = c(1, 4)
)
)
)
# Prepare plotting data
dims <- paste0(Key(object = object[[reduction]]), dims)
plot.data <- FetchData(object = object, vars = c(dims, feature), layer = slot)
# Shiny server
server <- function(input, output, session) {
plot.env <- reactiveValues(
data = plot.data,
dims = paste0(Key(object = object[[reduction]]), dims),
feature = feature,
palette = 'Seurat'
)
# Observe events
observeEvent(
eventExpr = input$done,
handlerExpr = stopApp(returnValue = plot.env$plot)
)
observe(x = {
assay <- input$assay
feature.use <- input$feature
features.assay <- sort(x = rownames(x = GetAssayData(
object = object,
layer = slot,
assay = assay
)))
feature.use <- ifelse(
test = feature.use %in% features.assay,
yes = feature.use,
no = features.assay[1]
)
reduc <- input$reduction
dims.reduc <- gsub(
pattern = Key(object = object[[reduc]]),
replacement = '',
x = colnames(x = object[[reduc]])
)
dims <- c(input$xdim, input$ydim)
for (i in seq_along(along.with = dims)) {
if (!dims[i] %in% dims.reduc) {
dims[i] <- dims.reduc[i]
}
}
updateSelectInput(
session = session,
inputId = 'xdim',
label = 'X dimension',
choices = dims.reduc,
selected = as.character(x = dims[1])
)
updateSelectInput(
session = session,
inputId = 'ydim',
label = 'Y dimension',
choices = dims.reduc,
selected = as.character(x = dims[2])
)
updateSelectInput(
session = session,
inputId = 'feature',
label = feature.label,
choices = features.assay,
selected = feature.use
)
})
observe(x = {
feature.use <- input$feature
feature.keyed <- paste0(Key(object = object[[input$assay]]), feature.use)
reduc <- input$reduction
dims <- c(input$xdim, input$ydim)
dims <- paste0(Key(object = object[[reduc]]), dims)
plot.data <- tryCatch(
expr = FetchData(
object = object,
vars = c(dims, feature.keyed),
layer = slot
),
warning = function(...) {
return(plot.env$data)
},
error = function(...) {
return(plot.env$data)
}
)
dims <- colnames(x = plot.data)[1:2]
colnames(x = plot.data) <- c(dims, feature.use)
plot.env$data <- plot.data
plot.env$feature <- feature.use
plot.env$dims <- dims
})
observe(x = {
plot.env$palette <- input$palette
})
# Create the plot
output$plot <- renderPlot(expr = {
plot.env$plot <- SingleDimPlot(
data = plot.env$data,
dims = plot.env$dims,
col.by = plot.env$feature,
label = FALSE
) +
theme_cowplot() +
theme(plot.title = element_text(hjust = 0.5)) +
guides(color = NULL) +
scale_color_gradientn(
colors = FeaturePalettes[[plot.env$palette]],
guide = 'colorbar'
)
plot.env$plot
})
}
runGadget(app = ui, server = server)
}
#' Highlight Neighbors in DimPlot
#'
#' It will color the query cells and the neighbors of the query cells in the
#' DimPlot
#'
#' @inheritParams DimPlot
#' @param nn.idx the neighbor index of all cells
#' @param query.cells cells used to find their neighbors
#' @param show.all.cells Show all cells or only query and neighbor cells
#'
#' @inherit DimPlot return
#'
#' @export
#' @concept visualization
#'
NNPlot <- function(
object,
reduction,
nn.idx,
query.cells,
dims = 1:2,
label = FALSE,
label.size = 4,
repel = FALSE,
sizes.highlight = 2,
pt.size = 1,
cols.highlight = c("#377eb8", "#e41a1c"),
na.value = "#bdbdbd",
order = c("self", "neighbors", "other"),
show.all.cells = TRUE,
...
) {
if (inherits(x = nn.idx, what = 'Neighbor')) {
rownames(x = slot(object = nn.idx, name = 'nn.idx')) <- Cells(x = nn.idx)
nn.idx <- Indices(object = nn.idx)
}
if (length(x = query.cells) > 1) {
neighbor.cells <- apply(
X = nn.idx[query.cells, -1],
MARGIN = 2,
FUN = function(x) {
return(Cells(x = object)[x])
}
)
} else {
neighbor.cells <- Cells(x = object)[nn.idx[query.cells , -1]]
}
neighbor.cells <- as.vector(x = neighbor.cells)
neighbor.cells <- neighbor.cells[!is.na(x = neighbor.cells)]
object[["nn.col"]] <- "other"
object[["nn.col"]][neighbor.cells, ] <- "neighbors"
object[["nn.col"]][query.cells, ] <- "self"
object$nn.col <- factor(
x = object$nn.col,
levels = c("self", "neighbors", "other")
)
if (!show.all.cells) {
object <- subset(
x = object,
cells = Cells(x = object)[which(x = object[["nn.col"]] != "other")]
)
nn.cols <- c(rev(x = cols.highlight))
nn.pt.size <- sizes.highlight
} else {
highlight.info <- SetHighlight(
cells.highlight = c(query.cells, neighbor.cells),
cells.all = Cells(x = object),
sizes.highlight = sizes.highlight,
pt.size = pt.size,
cols.highlight = "red"
)
nn.cols <- c(na.value, rev(x = cols.highlight))
nn.pt.size <- highlight.info$size
}
NN.plot <- DimPlot(
object = object,
reduction = reduction,
dims = dims,
group.by = "nn.col",
cols = nn.cols,
label = label,
order = order,
pt.size = nn.pt.size ,
label.size = label.size,
repel = repel
)
return(NN.plot)
}
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Scatter plots
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#' Cell-cell scatter plot
#'
#' Creates a plot of scatter plot of features across two single cells. Pearson
#' correlation between the two cells is displayed above the plot.
#'
#' @inheritParams FeatureScatter
#' @inheritParams DimPlot
#' @param cell1 Cell 1 name
#' @param cell2 Cell 2 name
#' @param features Features to plot (default, all features)
#' @param highlight Features to highlight
#' @return A ggplot object
#'
#' @export
#' @concept visualization
#'
#' @aliases CellPlot
#'
#' @examples
#' data("pbmc_small")
#' CellScatter(object = pbmc_small, cell1 = 'ATAGGAGAAACAGA', cell2 = 'CATCAGGATGCACA')
#'
CellScatter <- function(
object,
cell1,
cell2,
features = NULL,
highlight = NULL,
cols = NULL,
pt.size = 1,
smooth = FALSE,
raster = NULL,
raster.dpi = c(512, 512)
) {
features <- features %||% rownames(x = object)
data <- FetchData(
object = object,
vars = features,
cells = c(cell1, cell2)
)
data <- as.data.frame(x = t(x = data))
plot <- SingleCorPlot(
data = data,
cols = cols,
pt.size = pt.size,
rows.highlight = highlight,
smooth = smooth,
raster = raster,
raster.dpi = raster.dpi
)
return(plot)
}
#' Scatter plot of single cell data
#'
#' Creates a scatter plot of two features (typically feature expression), across a
#' set of single cells. Cells are colored by their identity class. Pearson
#' correlation between the two features is displayed above the plot.
#'
#' @param object Seurat object
#' @param feature1 First feature to plot. Typically feature expression but can also
#' be metrics, PC scores, etc. - anything that can be retreived with FetchData
#' @param feature2 Second feature to plot.
#' @param cells Cells to include on the scatter plot.
#' @param shuffle Whether to randomly shuffle the order of points. This can be
#' useful for crowded plots if points of interest are being buried. (default is FALSE)
#' @param seed Sets the seed if randomly shuffling the order of points.
#' @param group.by Name of one or more metadata columns to group (color) cells by
#' (for example, orig.ident); pass 'ident' to group by identity class
#' @param cols Colors to use for identity class plotting.
#' @param pt.size Size of the points on the plot
#' @param shape.by Ignored for now
#' @param split.by A factor in object metadata to split the feature plot by, pass 'ident'
#' to split by cell identity
#' @param span Spline span in loess function call, if \code{NULL}, no spline added
#' @param smooth Smooth the graph (similar to smoothScatter)
#' @param slot Slot to pull data from, should be one of 'counts', 'data', or 'scale.data'
#' @param combine Combine plots into a single \code{\link[patchwork]{patchwork}ed}
#' @param plot.cor Display correlation in plot title
#' @param ncol Number of columns if plotting multiple plots
#' @param raster Convert points to raster format, default is \code{NULL}
#' which will automatically use raster if the number of points plotted is greater than
#' 100,000
#' @param raster.dpi Pixel resolution for rasterized plots, passed to geom_scattermore().
#' Default is c(512, 512).
#' @param jitter Jitter for easier visualization of crowded points (default is FALSE)
#' @param log Plot features on the log scale (default is FALSE)
#'
#' @return A ggplot object
#'
#' @importFrom ggplot2 geom_smooth facet_wrap vars sym labs
#' @importFrom patchwork wrap_plots
#'
#' @export
#' @concept visualization
#'
#' @aliases GenePlot
#'
#' @examples
#' data("pbmc_small")
#' FeatureScatter(object = pbmc_small, feature1 = 'CD9', feature2 = 'CD3E')
#'
FeatureScatter <- function(
object,
feature1,
feature2,
cells = NULL,
shuffle = FALSE,
seed = 1,
group.by = NULL,
split.by = NULL,
cols = NULL,
pt.size = 1,
shape.by = NULL,
span = NULL,
smooth = FALSE,
combine = TRUE,
slot = 'data',
plot.cor = TRUE,
ncol = NULL,
raster = NULL,
raster.dpi = c(512, 512),
jitter = FALSE,
log = FALSE
) {
cells <- cells %||% colnames(x = object)
if (isTRUE(x = shuffle)) {
set.seed(seed = seed)
cells <- sample(x = cells)
}
group.by <- group.by %||% 'ident'
data <- FetchData(
object = object,
vars = c(feature1, feature2, group.by),
cells = cells,
layer = slot
)
if (!grepl(pattern = feature1, x = names(x = data)[1], fixed = TRUE)) {
abort(message = paste("Feature 1", sQuote(x = feature1), "not found"))
}
if (!grepl(pattern = feature2, x = names(x = data)[2], fixed = TRUE)) {
abort(message = paste("Feature 2", sQuote(x = feature2), "not found"))
}
feature1 <- names(x = data)[1]
feature2 <- names(x = data)[2]
group.by <- intersect(x = group.by, y = names(x = data)[3:ncol(x = data)])
for (group in group.by) {
if (!is.factor(x = data[, group])) {
data[, group] <- factor(x = data[, group])
}
}
if (!is.null(x = split.by)) {
split <- FetchData(object = object, vars = split.by, clean=TRUE)[split.by]
data <- data[rownames(split),]
data[, split.by] <- split
}
plots <- lapply(
X = group.by,
FUN = function(x) {
plot <- SingleCorPlot(
data = data[,c(feature1, feature2, split.by)],
col.by = data[, x],
cols = cols,
pt.size = pt.size,
smooth = smooth,
legend.title = 'Identity',
span = span,
plot.cor = plot.cor,
raster = raster,
raster.dpi = raster.dpi,
jitter = jitter
)
if (!is.null(x = split.by)) {
plot <- plot + FacetTheme() +
facet_wrap(
facets = vars(!!sym(x = split.by)),
ncol = if (length(x = group.by) > 1 || is.null(x = ncol)) {
length(x = unique(x = data[, split.by]))
} else {
ncol
}
)
}
if (log) {
plot <- plot + scale_x_log10() + scale_y_log10()
}
plot
}
)
if (isTRUE(x = length(x = plots) == 1)) {
return(plots[[1]])
}
if (isTRUE(x = combine)) {
plots <- wrap_plots(plots, ncol = length(x = group.by))
}
return(plots)
}
#' View variable features
#'
#' @inheritParams FeatureScatter
#' @inheritParams SeuratObject::HVFInfo
#' @param cols Colors to specify non-variable/variable status
#' @param assay Assay to pull variable features from
#' @param log Plot the x-axis in log scale
#' @param raster Convert points to raster format, default is \code{NULL}
#' which will automatically use raster if the number of points plotted is greater than
#' 100,000
#'
#' @return A ggplot object
#'
#' @importFrom ggplot2 labs scale_color_manual scale_x_log10
#' @export
#' @concept visualization
#'
#' @aliases VariableGenePlot MeanVarPlot
#'
#' @seealso \code{\link{FindVariableFeatures}}
#'
#' @examples
#' data("pbmc_small")
#' VariableFeaturePlot(object = pbmc_small)
#'
VariableFeaturePlot <- function(
object,
cols = c('black', 'red'),
pt.size = 1,
log = NULL,
selection.method = NULL,
assay = NULL,
raster = NULL,
raster.dpi = c(512, 512)
) {
if (length(x = cols) != 2) {
stop("'cols' must be of length 2")
}
hvf.info <- HVFInfo(
object = object,
assay = assay,
method = selection.method,
status = TRUE
)
status.col <- colnames(hvf.info)[grepl("variable", colnames(hvf.info))][[1]]
var.status <- c('no', 'yes')[unlist(hvf.info[[status.col]]) + 1]
if (colnames(x = hvf.info)[3] == 'dispersion.scaled') {
hvf.info <- hvf.info[, c(1, 2)]
} else if (colnames(x = hvf.info)[3] == 'variance.expected') {
hvf.info <- hvf.info[, c(1, 4)]
} else {
hvf.info <- hvf.info[, c(1, 3)]
}
axis.labels <- switch(
EXPR = colnames(x = hvf.info)[2],
'variance.standardized' = c('Average Expression', 'Standardized Variance'),
'dispersion' = c('Average Expression', 'Dispersion'),
'residual_variance' = c('Geometric Mean of Expression', 'Residual Variance')
)
log <- log %||% (any(c('variance.standardized', 'residual_variance') %in% colnames(x = hvf.info)))
# var.features <- VariableFeatures(object = object, assay = assay)
# var.status <- ifelse(
# test = rownames(x = hvf.info) %in% var.features,
# yes = 'yes',
# no = 'no'
# )
plot <- SingleCorPlot(
data = hvf.info,
col.by = var.status,
pt.size = pt.size,
raster = raster,
raster.dpi = raster.dpi
)
if (length(x = unique(x = var.status)) == 1) {
switch(
EXPR = var.status[1],
'yes' = {
cols <- cols[2]
labels.legend <- 'Variable'
},
'no' = {
cols <- cols[1]
labels.legend <- 'Non-variable'
}
)
} else {
labels.legend <- c('Non-variable', 'Variable')
}
plot <- plot +
labs(title = NULL, x = axis.labels[1], y = axis.labels[2]) +
scale_color_manual(
labels = paste(labels.legend, 'count:', table(var.status)),
values = cols
)
if (log) {
plot <- plot + scale_x_log10()
}
return(plot)
}
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Polygon Plots
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#' Polygon DimPlot
#'
#' Plot cells as polygons, rather than single points. Color cells by identity, or a categorical variable
#' in metadata
#'
#' @inheritParams PolyFeaturePlot
#' @param group.by A grouping variable present in the metadata. Default is to use the groupings present
#' in the current cell identities (\code{Idents(object = object)})
#'
#' @return Returns a ggplot object
#'
#' @export
#' @concept visualization
#'
PolyDimPlot <- function(
object,
group.by = NULL,
cells = NULL,
poly.data = 'spatial',
flip.coords = FALSE
) {
polygons <- Misc(object = object, slot = poly.data)
if (is.null(x = polygons)) {
stop("Could not find polygon data in misc slot")
}
group.by <- group.by %||% 'ident'
group.data <- FetchData(
object = object,
vars = group.by,
cells = cells
)
group.data$cell <- rownames(x = group.data)
data <- merge(x = polygons, y = group.data, by = 'cell')
if (flip.coords) {
coord.x <- data$x
data$x <- data$y
data$y <- coord.x
}
plot <- SinglePolyPlot(data = data, group.by = group.by)
return(plot)
}
#' Polygon FeaturePlot
#'
#' Plot cells as polygons, rather than single points. Color cells by any value
#' accessible by \code{\link{FetchData}}.
#'
#' @inheritParams FeaturePlot
#' @param poly.data Name of the polygon dataframe in the misc slot
#' @param ncol Number of columns to split the plot into
#' @param common.scale ...
#' @param flip.coords Flip x and y coordinates
#'
#' @return Returns a ggplot object
#'
#' @importFrom ggplot2 scale_fill_viridis_c facet_wrap
#'
#' @export
#' @concept visualization
#' @concept spatial
#'
PolyFeaturePlot <- function(
object,
features,
cells = NULL,
poly.data = 'spatial',
ncol = ceiling(x = length(x = features) / 2),
min.cutoff = 0,
max.cutoff = NA,
common.scale = TRUE,
flip.coords = FALSE
) {
polygons <- Misc(object = object, slot = poly.data)
if (is.null(x = polygons)) {
stop("Could not find polygon data in misc slot")
}
assay.data <- FetchData(
object = object,
vars = features,
cells = cells
)
features <- colnames(x = assay.data)
cells <- rownames(x = assay.data)
min.cutoff <- mapply(
FUN = function(cutoff, feature) {
return(ifelse(
test = is.na(x = cutoff),
yes = min(assay.data[, feature]),
no = cutoff
))
},
cutoff = min.cutoff,
feature = features
)
max.cutoff <- mapply(
FUN = function(cutoff, feature) {
return(ifelse(
test = is.na(x = cutoff),
yes = max(assay.data[, feature]),
no = cutoff
))
},
cutoff = max.cutoff,
feature = features
)
check.lengths <- unique(x = vapply(
X = list(features, min.cutoff, max.cutoff),
FUN = length,
FUN.VALUE = numeric(length = 1)
))
if (length(x = check.lengths) != 1) {
stop("There must be the same number of minimum and maximum cuttoffs as there are features")
}
assay.data <- mapply(
FUN = function(feature, min, max) {
return(ScaleColumn(vec = assay.data[, feature], cutoffs = c(min, max)))
},
feature = features,
min = min.cutoff,
max = max.cutoff
)
if (common.scale) {
assay.data <- apply(
X = assay.data,
MARGIN = 2,
FUN = function(x) {
return(x - min(x))
}
)
assay.data <- t(
x = t(x = assay.data) / apply(X = assay.data, MARGIN = 2, FUN = max)
)
}
assay.data <- as.data.frame(x = assay.data)
assay.data <- data.frame(
cell = as.vector(x = replicate(n = length(x = features), expr = cells)),
feature = as.vector(x = t(x = replicate(n = length(x = cells), expr = features))),
expression = unlist(x = assay.data, use.names = FALSE)
)
data <- merge(x = polygons, y = assay.data, by = 'cell')
data$feature <- factor(x = data$feature, levels = features)
if (flip.coords) {
coord.x <- data$x
data$x <- data$y
data$y <- coord.x
}
plot <- SinglePolyPlot(data = data, group.by = 'expression', font_size = 8) +
scale_fill_viridis_c() +
facet_wrap(facets = 'feature', ncol = ncol)
return(plot)
}
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Spatial Plots
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#' Spatial Cluster Plots
#'
#' Visualize clusters or other categorical groupings in a spatial context
#'
#' @inheritParams DimPlot
#' @inheritParams SingleImagePlot
#' @param object A \code{\link[SeuratObject]{Seurat}} object
#' @param fov Name of FOV to plot
#' @param boundaries A vector of segmentation boundaries per image to plot;
#' can be a character vector, a named character vector, or a named list.
#' Names should be the names of FOVs and values should be the names of
#' segmentation boundaries
#' @param molecules A vector of molecules to plot
#' @param nmols Max number of each molecule specified in `molecules` to plot
#' @param dark.background Set plot background to black
#' @param crop Crop the plots to area with cells only
#' @param overlap Overlay boundaries from a single image to create a single
#' plot; if \code{TRUE}, then boundaries are stacked in the order they're
#' given (first is lowest)
#' @param axes Keep axes and panel background
#' @param combine Combine plots into a single
#' \code{patchwork} ggplot object.If \code{FALSE},
#' return a list of ggplot objects
#' @param coord.fixed Plot cartesian coordinates with fixed aspect ratio
#' @param flip_xy Flag to flip X and Y axes. Default is FALSE.
#'
#' @return If \code{combine = TRUE}, a \code{patchwork}
#' ggplot object; otherwise, a list of ggplot objects
#'
#' @importFrom rlang !! is_na sym
#' @importFrom patchwork wrap_plots
#' @importFrom ggplot2 element_blank facet_wrap vars
#' @importFrom SeuratObject DefaultFOV Cells
#' DefaultBoundary FetchData Images Overlay
#'
#' @export
#' @concept visualization
#' @concept spatial
#'
ImageDimPlot <- function(
object,
fov = NULL,
boundaries = NULL,
group.by = NULL,
split.by = NULL,
cols = NULL,
shuffle.cols = FALSE,
size = 0.5,
molecules = NULL,
mols.size = 0.1,
mols.cols = NULL,
mols.alpha = 1.0,
nmols = 1000,
alpha = 1.0,
border.color = 'white',
border.size = NULL,
na.value = 'grey50',
dark.background = TRUE,
crop = FALSE,
cells = NULL,
overlap = FALSE,
axes = FALSE,
combine = TRUE,
coord.fixed = TRUE,
flip_xy = TRUE
) {
cells <- cells %||% Cells(x = object)
# Determine FOV to use
fov <- fov %||% DefaultFOV(object = object)
fov <- Filter(
f = function(x) {
return(
x %in% Images(object = object) &&
inherits(x = object[[x]], what = 'FOV')
)
},
x = fov
)
if (!length(x = fov)) {
stop("No compatible spatial coordinates present")
}
# Identify boundaries to use
boundaries <- boundaries %||% sapply(
X = fov,
FUN = function(x) {
return(DefaultBoundary(object = object[[x]]))
},
simplify = FALSE,
USE.NAMES = TRUE
)
boundaries <- .BoundariesByImage(
object = object,
fov = fov,
boundaries = boundaries
)
fov <- names(x = boundaries)
overlap <- rep_len(x = overlap, length.out = length(x = fov))
crop <- rep_len(x = crop, length.out = length(x = fov))
names(x = crop) <- fov
# Prepare plotting data
group.by <- boundaries %!NA% group.by %||% 'ident'
vars <- c(group.by, split.by)
md <- if (!is_na(x = vars)) {
FetchData(
object = object,
vars = vars[!is.na(x = vars)],
cells = cells
)
} else {
NULL
}
pnames <- unlist(x = lapply(
X = seq_along(along.with = fov),
FUN = function(i) {
return(if (isTRUE(x = overlap[i])) {
fov[i]
} else {
paste(fov[i], boundaries[[i]], sep = '_')
})
}
))
pdata <- vector(mode = 'list', length = length(x = pnames))
names(x = pdata) <- pnames
for (i in names(x = pdata)) {
ul <- unlist(x = strsplit(x = i, split = '_'))
img <- paste(ul[1:length(ul)-1], collapse = '_')
# Apply overlap
lyr <- ul[length(ul)]
if (is.na(x = lyr)) {
lyr <- boundaries[[img]]
}
# TODO: Apply crop
pdata[[i]] <- lapply(
X = lyr,
FUN = function(l) {
if (l == 'NA') {
return(NA)
}
df <- fortify(model = object[[img]][[l]])
df <- df[df$cell %in% cells, , drop = FALSE]
if (!is.null(x = md)) {
df <- merge(x = df, y = md, by.x = 'cell', by.y = 0, all.x = TRUE)
}
df$cell <- paste(l, df$cell, sep = '_')
df$boundary <- l
return(df)
}
)
pdata[[i]] <- if (!is_na(x = pdata[[i]])) {
do.call(what = 'rbind', args = pdata[[i]])
} else {
unlist(x = pdata[[i]])
}
}
# Fetch molecule information
if (!is.null(x = molecules)) {
molecules <- .MolsByFOV(
object = object,
fov = fov,
molecules = molecules
)
mdata <- vector(mode = 'list', length = length(x = fov))
names(x = mdata) <- fov
for (img in names(x = mdata)) {
idata <- object[[img]]
if (!img %in% names(x = molecules)) {
mdata[[img]] <- NULL
next
}
if (isTRUE(x = crop[img])) {
idata <- Overlay(x = idata, y = idata)
}
imols <- gsub(
pattern = paste0('^', Key(object = idata)),
replacement = '',
x = molecules[[img]]
)
mdata[[img]] <- FetchData(
object = idata,
vars = imols,
nmols = nmols
)
}
} else {
mdata <- NULL
}
# Build the plots
plots <- vector(
mode = 'list',
length = length(x = pdata) * ifelse(
test = length(x = group.by),
yes = length(x = group.by),
no = 1L
)
)
idx <- 1L
for (group in group.by) {
for (i in seq_along(along.with = pdata)) {
img <- unlist(x = strsplit(x = names(x = pdata)[i], split = '_'))[1L]
p <- SingleImagePlot(
data = pdata[[i]],
col.by = pdata[[i]] %!NA% group,
molecules = mdata[[img]],
cols = cols,
shuffle.cols = shuffle.cols,
size = size,
alpha = alpha,
mols.size = mols.size,
mols.cols = mols.cols,
mols.alpha = mols.alpha,
border.color = border.color,
border.size = border.size,
na.value = na.value,
dark.background = dark.background
)
if (!is.null(x = split.by)) {
p <- p + facet_wrap(
facets = vars(!!sym(x = split.by))
)
}
if (!isTRUE(x = axes)) {
p <- p + NoAxes(panel.background = element_blank())
}
if (!anyDuplicated(x = pdata[[i]]$cell)) {
p <- p + guides(fill = guide_legend(override.aes = list(size=4L, alpha=1)))
}
if (isTRUE(coord.fixed)) {
p <- p + coord_fixed()
}
if(!isTRUE(flip_xy) && isTRUE(coord.fixed)){
xy_ratio = (max(pdata[[i]]$x) - min(pdata[[i]]$x)) / (max(pdata[[i]]$y) - min(pdata[[i]]$y))
p = p + coord_flip() + theme(aspect.ratio = 1/xy_ratio)
}
plots[[idx]] <- p
idx <- idx + 1L
}
}
if (isTRUE(x = combine)) {
plots <- wrap_plots(plots)
}
return(plots)
}
#' Spatial Feature Plots
#'
#' Visualize expression in a spatial context
#'
#' @inheritParams FeaturePlot
#' @inheritParams ImageDimPlot
#' @param scale Set color scaling across multiple plots; choose from:
#' \itemize{
#' \item \dQuote{\code{feature}}: Plots per-feature are scaled across splits
#' \item \dQuote{\code{all}}: Plots per-feature are scaled across all features
#' \item \dQuote{\code{none}}: Plots are not scaled; \strong{note}: setting
#' \code{scale} to \dQuote{\code{none}} will result in color scales that are
#' \emph{not} comparable between plots
#' }
#' Ignored if \code{blend = TRUE}
#'
#' @inherit ImageDimPlot return
#'
#' @importFrom patchwork wrap_plots
#' @importFrom cowplot theme_cowplot
#' @importFrom ggplot2 dup_axis element_blank element_text facet_wrap guides
#' labs margin vars scale_y_continuous theme
#' @importFrom SeuratObject DefaultFOV Cells DefaultBoundary
#' FetchData Images Overlay
#'
#' @export
#' @concept visualization
#' @concept spatial
#'
ImageFeaturePlot <- function(
object,
features,
fov = NULL,
boundaries = NULL,
cols = if (isTRUE(x = blend)) {
c("lightgrey", "#ff0000", "#00ff00")
} else {
c("lightgrey", "firebrick1")
},
size = 0.5,
min.cutoff = NA,
max.cutoff = NA,
split.by = NULL,
molecules = NULL,
mols.size = 0.1,
mols.cols = NULL,
nmols = 1000,
alpha = 1.0,
border.color = 'white',
border.size = NULL,
dark.background = TRUE,
blend = FALSE,
blend.threshold = 0.5,
crop = FALSE,
cells = NULL,
scale = c('feature', 'all', 'none'),
overlap = FALSE,
axes = FALSE,
combine = TRUE,
coord.fixed = TRUE
) {
cells <- cells %||% Cells(x = object)
scale <- scale[[1L]]
scale <- match.arg(arg = scale)
# Set a theme to remove right-hand Y axis lines
# Also sets right-hand Y axis text label formatting
no.right <- theme(
axis.line.y.right = element_blank(),
axis.ticks.y.right = element_blank(),
axis.text.y.right = element_blank(),
axis.title.y.right = element_text(
face = "bold",
size = 14,
margin = margin(r = 7)
)
)
# Determine fov to use
fov <- fov %||% DefaultFOV(object = object)
fov <- Filter(
f = function(x) {
return(
x %in% Images(object = object) &&
inherits(x = object[[x]], what = 'FOV')
)
},
x = fov
)
if (!length(x = fov)) {
stop("No compatible spatial coordinates present")
}
# Identify boundaries to use
boundaries <- boundaries %||% sapply(
X = fov,
FUN = function(x) {
return(DefaultBoundary(object = object[[x]]))
},
simplify = FALSE,
USE.NAMES = TRUE
)
boundaries <- .BoundariesByImage(
object = object,
fov = fov,
boundaries = boundaries
)
fov <- names(x = boundaries)
# Check overlaps/crops
if (isTRUE(x = blend) || !is.null(x = split.by)) {
type <- ifelse(test = isTRUE(x = 'blend'), yes = 'Blended', no = 'Split')
if (length(x = fov) != 1L) {
fov <- fov[1L]
warning(
type,
' image feature plots can only be done on a single image, using "',
fov,
'"',
call. = FALSE,
immediate. = TRUE
)
}
if (any(!overlap) && length(x = boundaries[[fov]]) > 1L) {
warning(
type,
" image feature plots require overlapped segmentations",
call. = FALSE,
immediate. = TRUE
)
}
overlap <- TRUE
}
overlap <- rep_len(x = overlap, length.out = length(x = fov))
crop <- rep_len(x = crop, length.out = length(x = fov))
names(x = crop) <- names(x = overlap) <- fov
# Checks for blending
if (isTRUE(x = blend)) {
if (length(x = features) != 2L) {
stop("Blended feature plots only works with two features")
}
default.colors <- eval(expr = formals(fun = ImageFeaturePlot)$cols)
cols <- switch(
EXPR = as.character(x = length(x = cols)),
'0' = {
warning("No colors provided, using default colors", immediate. = TRUE)
default.colors
},
'1' = {
warning(
"Only one color provided, assuming specified is double-negative and augmenting with default colors",
immediate. = TRUE
)
c(cols, default.colors[2:3])
},
'2' = {
warning(
"Only two colors provided, assuming specified are for features and augmenting with '",
default.colors[1],
"' for double-negatives",
immediate. = TRUE
)
c(default.colors[1], cols)
},
'3' = cols,
{
warning(
"More than three colors provided, using only first three",
immediate. = TRUE
)
cols[1:3]
}
)
}
# Get feature, splitting data
md <- FetchData(
object = object,
vars = c(features, split.by[1L]),
cells = cells
)
split.by <- intersect(x = split.by, y = colnames(x = md))
if (!length(x = split.by)) {
split.by <- NULL
}
imax <- ifelse(
test = is.null(x = split.by),
yes = ncol(x = md),
no = ncol(x = md) - length(x = split.by)
)
features <- colnames(x = md)[1:imax]
# Determine cutoffs
min.cutoff <- mapply(
FUN = function(cutoff, feature) {
return(ifelse(
test = is.na(x = cutoff),
yes = min(md[[feature]]),
no = cutoff
))
},
cutoff = min.cutoff,
feature = features
)
max.cutoff <- mapply(
FUN = function(cutoff, feature) {
return(ifelse(
test = is.na(x = cutoff),
yes = max(md[[feature]]),
no = cutoff
))
},
cutoff = max.cutoff,
feature = features
)
check.lengths <- unique(x = vapply(
X = list(features, min.cutoff, max.cutoff),
FUN = length,
FUN.VALUE = numeric(length = 1)
))
if (length(x = check.lengths) != 1) {
stop("There must be the same number of minimum and maximum cuttoffs as there are features")
}
brewer.gran <- ifelse(
test = length(x = cols) == 1,
yes = brewer.pal.info[cols, ]$maxcolors,
no = length(x = cols)
)
# Apply cutoffs
for (i in seq_along(along.with = features)) {
f <- features[[i]]
data.feature <- md[[f]]
min.use <- SetQuantile(cutoff = min.cutoff[i], data = data.feature)
max.use <- SetQuantile(cutoff = max.cutoff[i], data = data.feature)
data.feature[data.feature < min.use] <- min.use
data.feature[data.feature > max.use] <- max.use
if (brewer.gran != 2) {
data.feature <- if (all(data.feature == 0)) {
rep_len(x = 0, length.out = length(x = data.feature))
} else {
as.numeric(x = as.factor(x = cut(
x = as.numeric(x = data.feature),
breaks = brewer.gran
)))
}
}
md[[f]] <- data.feature
}
# Figure out splits
if (is.null(x = split.by)) {
split.by <- RandomName()
md[[split.by]] <- factor(x = split.by)
}
if (!is.factor(x = md[[split.by]])) {
md[[split.by]] <- factor(x = md[[split.by]])
}
# Apply blends
if (isTRUE(x = blend)) {
md <- lapply(
X = levels(x = md[[split.by]]),
FUN = function(x) {
df <- md[as.character(x = md[[split.by]]) == x, , drop = FALSE]
no.expression <- features[colMeans(x = df[, features]) == 0]
if (length(x = no.expression)) {
stop(
"The following features have no value: ",
paste(no.expression, collapse = ', ')
)
}
return(cbind(
df[, split.by, drop = FALSE],
BlendExpression(data = df[, features])
))
}
)
md <- do.call(what = 'rbind', args = md)
features <- setdiff(x = colnames(x = md), y = split.by)
}
# Prepare plotting data
pnames <- unlist(x = lapply(
X = seq_along(along.with = fov),
FUN = function(i) {
return(if (isTRUE(x = overlap[i])) {
fov[i]
} else {
paste(fov[i], boundaries[[i]], sep = '_')
})
}
))
pdata <- vector(mode = 'list', length = length(x = pnames))
names(x = pdata) <- pnames
for (i in names(x = pdata)) {
ul <- unlist(x = strsplit(x = i, split = '_'))
# img <- paste(ul[1:length(ul)-1], collapse = '_')
# Apply overlap
# lyr <- ul[length(ul)]
if(length(ul) > 1) {
img <- paste(ul[1:length(ul)-1], collapse = '_')
lyr <- ul[length(ul)]
} else if (length(ul) == 1) {
img <- ul[1]
lyr <- "centroids"
} else {
stop("the length of ul is 0. please check.")
}
if (is.na(x = lyr)) {
lyr <- boundaries[[img]]
}
pdata[[i]] <- lapply(
X = lyr,
FUN = function(l) {
df <- fortify(model = object[[img]][[l]])
df <- df[df$cell %in% cells, , drop = FALSE]
if (!is.null(x = md)) {
df <- merge(x = df, y = md, by.x = 'cell', by.y = 0, all.x = TRUE)
}
df$cell <- paste(l, df$cell, sep = '_')
df$boundary <- l
return(df)
}
)
pdata[[i]] <- if (!is_na(x = pdata[[i]])) {
do.call(what = 'rbind', args = pdata[[i]])
} else {
unlist(x = pdata[[i]])
}
}
# Fetch molecule information
if (!is.null(x = molecules)) {
molecules <- .MolsByFOV(
object = object,
fov = fov,
molecules = molecules
)
mdata <- vector(mode = 'list', length = length(x = fov))
names(x = mdata) <- fov
for (img in names(x = mdata)) {
idata <- object[[img]]
if (!img %in% names(x = molecules)) {
mdata[[img]] <- NULL
next
}
if (isTRUE(x = crop[img])) {
idata <- Overlay(x = idata, y = idata)
}
imols <- gsub(
pattern = paste0('^', Key(object = idata)),
replacement = '',
x = molecules[[img]]
)
mdata[[img]] <- FetchData(
object = idata,
vars = imols,
nmols = nmols
)
}
} else {
mdata <- NULL
}
# Set blended colors
if (isTRUE(x = blend)) {
ncol <- 4
color.matrix <- BlendMatrix(
two.colors = cols[2:3],
col.threshold = blend.threshold,
negative.color = cols[1]
)
cols <- cols[2:3]
colors <- list(
color.matrix[, 1],
color.matrix[1, ],
as.vector(x = color.matrix)
)
blend.legend <- BlendMap(color.matrix = color.matrix)
}
limits <- switch(
EXPR = scale,
'all' = range(unlist(x = md[, features])),
NULL
)
# Build the plots
plots <- vector(
mode = 'list',
length = length(x = levels(x = md[[split.by]]))
)
names(x = plots) <- levels(x = md[[split.by]])
for (i in seq_along(along.with = levels(x = md[[split.by]]))) {
ident <- levels(x = md[[split.by]])[i]
plots[[ident]] <- vector(mode = 'list', length = length(x = pdata))
names(x = plots[[ident]]) <- names(x = pdata)
if (isTRUE(x = blend)) {
blend.key <- suppressMessages(
expr = blend.legend +
scale_y_continuous(
sec.axis = dup_axis(name = ifelse(
test = length(x = levels(x = md[[split.by]])) > 1,
yes = ident,
no = ''
)),
expand = c(0, 0)
) +
labs(
x = features[1L],
y = features[2L],
title = if (i == 1L) {
paste('Color threshold:', blend.threshold)
} else {
NULL
}
) +
no.right
)
}
for (j in seq_along(along.with = pdata)) {
key <- names(x = pdata)[j]
img <- unlist(x = strsplit(x = key, split = '_'))[1L]
plots[[ident]][[key]] <- vector(
mode = 'list',
length = length(x = features) + ifelse(
test = isTRUE(x = blend),
yes = 1L,
no = 0L
)
)
data.plot <- pdata[[j]][as.character(x = pdata[[j]][[split.by]]) == ident, , drop = FALSE]
for (y in seq_along(along.with = features)) {
feature <- features[y]
# Get blended colors
cols.use <- if (isTRUE(x = blend)) {
cc <- as.numeric(x = as.character(x = data.plot[, feature])) + 1
colors[[y]][sort(unique(x = cc))]
} else {
NULL
}
colnames(data.plot) <- gsub("-", "_", colnames(data.plot))
p <- SingleImagePlot(
data = data.plot,
col.by = gsub("-", "_", feature),
size = size,
col.factor = blend,
cols = cols.use,
molecules = mdata[[img]],
mols.size = mols.size,
mols.cols = mols.cols,
alpha = alpha,
border.color = border.color,
border.size = border.size,
dark.background = dark.background
) +
CenterTitle() + labs(fill=feature)
# Remove fill guides for blended plots
if (isTRUE(x = blend)) {
p <- p + guides(fill = 'none')
}
if (isTRUE(coord.fixed)) {
p <- p + coord_fixed()
}
# Remove axes
if (!isTRUE(x = axes)) {
p <- p + NoAxes(panel.background = element_blank())
} else if (isTRUE(x = blend) || length(x = levels(x = md[[split.by]])) > 1L) {
if (y != 1L) {
p <- p + theme(
axis.line.y = element_blank(),
axis.ticks.y = element_blank(),
axis.text.y = element_blank(),
axis.title.y.left = element_blank()
)
}
if (i != length(x = levels(x = md[[split.by]]))) {
p <- p + theme(
axis.line.x = element_blank(),
axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
axis.title.x = element_blank()
)
}
}
# Add colors for unblended plots
if (!isTRUE(x = blend)) {
if (length(x = cols) == 1L) {
p <- p + scale_fill_brewer(palette = cols)
} else {
cols.grad <- cols
fexp <- data.plot[data.plot[[split.by]] == ident, feature, drop = TRUE]
fexp <- unique(x = fexp)
if (length(x = fexp) == 1L) {
warning(
"All cells have the same value (",
fexp,
") of ",
feature,
call. = FALSE,
immediate. = TRUE
)
if (fexp == 0) {
cols.grad <- cols.grad[1L]
}
}
# Check if we're scaling the colorbar across splits
if (scale == 'feature') {
limits <- range(pdata[[j]][[feature]])
}
p <- p + ggplot2::scale_fill_gradientn(
colors = cols.grad,
guide = 'colorbar',
limits = limits
)
}
}
# Add some labels
p <- p + if (i == 1L) {
ggplot2::labs(title = feature)
} else {
ggplot2::labs(title = NULL)
}
plots[[ident]][[key]][[y]] <- p
}
if (isTRUE(x = blend)) {
plots[[ident]][[key]][[length(x = plots[[ident]][[key]])]] <- blend.key
} else if (length(x = levels(x = md[[split.by]])) > 1L) {
plots[[ident]][[key]][[y]] <- suppressMessages(
expr = plots[[ident]][[key]][[y]] +
scale_y_continuous(sec.axis = dup_axis(name = ident)) +
no.right
)
}
}
plots[[ident]] <- unlist(
x = plots[[ident]],
recursive = FALSE,
use.names = FALSE
)
}
plots <- unlist(x = plots, recursive = FALSE, use.names = FALSE)
if (isTRUE(x = combine)) {
if (isTRUE(x = blend) || length(x = levels(x = md[[split.by]])) > 1L) {
plots <- wrap_plots(
plots,
ncol = ifelse(
test = isTRUE(x = blend),
yes = 4L,
no = length(x = features)
),
nrow = length(x = levels(x = md[[split.by]])),
guides = 'collect'
)
} else {
plots <- wrap_plots(plots)
}
}
return(plots)
}
#' Visualize spatial and clustering (dimensional reduction) data in a linked,
#' interactive framework
#'
#' @inheritParams SpatialPlot
#' @inheritParams FeaturePlot
#' @inheritParams DimPlot
#' @param feature Feature to visualize
#' @param image Name of the image to use in the plot
#'
#' @return Returns final plots. If \code{combine}, plots are stiched together
#' using \code{\link{CombinePlots}}; otherwise, returns a list of ggplot objects
#'
#' @rdname LinkedPlots
#' @name LinkedPlots
#'
#' @importFrom scales hue_pal
#' @importFrom patchwork wrap_plots
#' @importFrom ggplot2 scale_alpha_ordinal guides
#' @importFrom miniUI miniPage gadgetTitleBar miniTitleBarButton miniContentPanel
#' @importFrom shiny fillRow plotOutput brushOpts clickOpts hoverOpts
#' verbatimTextOutput reactiveValues observeEvent stopApp nearPoints
#' brushedPoints renderPlot renderPrint runGadget
#'
#' @aliases LinkedPlot LinkedDimPlot
#'
#' @export
#' @concept visualization
#' @concept spatial
#'
#' @examples
#' \dontrun{
#' LinkedDimPlot(seurat.object)
#' LinkedFeaturePlot(seurat.object, feature = 'Hpca')
#' }
#'
LinkedDimPlot <- function(
object,
dims = 1:2,
reduction = NULL,
image = NULL,
image.scale = "lowres",
group.by = NULL,
alpha = c(0.1, 1),
combine = TRUE
) {
# Setup gadget UI
ui <- miniPage(
gadgetTitleBar(
title = 'LinkedDimPlot',
left = miniTitleBarButton(inputId = 'reset', label = 'Reset')
),
miniContentPanel(
fillRow(
plotOutput(
outputId = 'spatialplot',
height = '100%',
# brush = brushOpts(id = 'brush', delay = 10, clip = TRUE, resetOnNew = FALSE),
click = clickOpts(id = 'spclick', clip = TRUE),
hover = hoverOpts(id = 'sphover', delay = 10, nullOutside = TRUE)
),
plotOutput(
outputId = 'dimplot',
height = '100%',
brush = brushOpts(id = 'brush', delay = 10, clip = TRUE, resetOnNew = FALSE),
click = clickOpts(id = 'dimclick', clip = TRUE),
hover = hoverOpts(id = 'dimhover', delay = 10, nullOutside = TRUE)
),
height = '97%'
),
verbatimTextOutput(outputId = 'info')
)
)
# Prepare plotting data
image <- image %||% DefaultImage(object = object)
cells.use <- Cells(x = object[[image]])
reduction <- reduction %||% DefaultDimReduc(object = object)
dims <- dims[1:2]
dims <- paste0(Key(object = object[[reduction]]), dims)
group.by <- group.by %||% 'ident'
group.data <- FetchData(
object = object,
vars = group.by,
cells = cells.use
)
coords <- GetTissueCoordinates(object = object[[image]], scale = image.scale)
embeddings <- Embeddings(object = object[[reduction]])[cells.use, dims]
plot.data <- cbind(coords, group.data, embeddings)
plot.data$selected_ <- FALSE
Idents(object = object) <- group.by
# Retrieve coordinates for tissue plot and dim plot separately
sp_x <- colnames(coords)[1]
sp_y <- colnames(coords)[2]
dp_x <- dims[1]
dp_y <- dims[2]
sp_y_min <- min(plot.data[[sp_y]])
sp_y_max <- max(plot.data[[sp_y]])
# Add tiny helper function to flip interactive coordinate points
flip_y <- function(pt) { if (!is.null(pt$y)) { pt$y <- sp_y_max - (pt$y - sp_y_min) }; pt }
# Setup the server
server <- function(input, output, session) {
click <- reactiveValues(pt = NULL, invert = FALSE)
plot.env <- reactiveValues(data = plot.data, alpha.by = NULL)
# Handle events
observeEvent(
eventExpr = input$done,
handlerExpr = {
plots <- list(plot.env$spatialplot, plot.env$dimplot)
if (combine) {
plots <- wrap_plots(plots, ncol = 2)
}
stopApp(returnValue = plots)
}
)
observeEvent(
eventExpr = input$reset,
handlerExpr = {
click$pt <- NULL
click$invert <- FALSE
plot.env$data <- plot.data
plot.env$alpha.by <- NULL
session$resetBrush(brushId = 'brush')
}
)
observeEvent(eventExpr = input$brush, handlerExpr = click$pt <- NULL)
observeEvent(
eventExpr = input$spclick,
handlerExpr = {
# Flip coordinates vertically for spatial plot to match tissue image
click$pt <- flip_y(input$spclick)
click$invert <- TRUE
}
)
observeEvent(
eventExpr = input$dimclick,
handlerExpr = {
click$pt <- input$dimclick
click$invert <- FALSE
}
)
observeEvent(
eventExpr = c(input$brush, input$spclick, input$dimclick),
handlerExpr = {
plot.env$data <- if (is.null(x = input$brush)) {
clicked <- nearPoints(
df = plot.data,
coordinfo = click$pt,
threshold = 10,
maxpoints = 1,
xvar = if (click$invert) sp_x else dp_x,
yvar = if (click$invert) sp_y else dp_y
)
if (nrow(x = clicked) == 1) {
cell.clicked <- rownames(x = clicked)
group.clicked <- plot.data[cell.clicked, group.by, drop = TRUE]
idx.group <- which(x = plot.data[[group.by]] == group.clicked)
plot.data[idx.group, 'selected_'] <- TRUE
plot.data
} else {
plot.data
}
} else if (input$brush$outputId == 'dimplot') {
brushedPoints(df = plot.data, brush = input$brush, allRows = TRUE, xvar = dp_x, yvar = dp_y)
} else if (input$brush$outputId == 'spatialplot') {
b <- input$brush
b$ymin <- sp_y_max - (b$ymin - sp_y_min)
b$ymax <- sp_y_max - (b$ymax - sp_y_min)
brushedPoints(df = plot.data, brush = b, allRows = TRUE, xvar = sp_x, yvar = sp_y)
}
plot.env$alpha.by <- if (any(plot.env$data$selected_)) {
'selected_'
} else {
NULL
}
}
)
# Set plots
output$spatialplot <- renderPlot(
expr = {
plot.env$spatialplot <- SingleSpatialPlot(
data = plot.env$data,
image = object[[image]],
col.by = group.by,
pt.size.factor = 1.6,
crop = TRUE,
alpha.by = plot.env$alpha.by
) + scale_alpha_ordinal(range = alpha) + NoLegend()
plot.env$spatialplot
}
)
output$dimplot <- renderPlot(
expr = {
plot.env$dimplot <- SingleDimPlot(
data = plot.env$data,
dims = dims,
col.by = group.by,
alpha.by = plot.env$alpha.by
) + scale_alpha_ordinal(range = alpha) + guides(alpha = "none")
plot.env$dimplot
}
)
# Add hover text
output$info <- renderPrint(
expr = {
cell.hover <- rownames(x = nearPoints(
df = plot.data,
coordinfo = if (is.null(input[['sphover']])) {
input$dimhover
} else {
flip_y(input$sphover)
},
threshold = 10,
maxpoints = 1,
xvar = if (is.null(input$sphover)) dp_x else sp_x,
yvar = if (is.null(input$sphover)) dp_y else sp_y
))
# if (length(x = cell.hover) == 1) {
# palette <- hue_pal()(n = length(x = levels(x = object)))
# group <- plot.data[cell.hover, group.by, drop = TRUE]
# background <- palette[which(x = levels(x = object) == group)]
# text <- unname(obj = BGTextColor(background = background))
# style <- paste0(
# paste(
# paste('background-color:', background),
# paste('color:', text),
# sep = '; '
# ),
# ';'
# )
# info <- paste(cell.hover, paste('Group:', group), sep = '<br />')
# } else {
# style <- 'background-color: white; color: black'
# info <- NULL
# }
# HTML(text = paste0("<div style='", style, "'>", info, "</div>"))
# p(HTML(info), style = style)
# paste0('<div style="', style, '">', info, '</div>')
# TODO: Get newlines, extra information, and background color working
if (length(x = cell.hover) == 1) {
paste(cell.hover, paste('Group:', plot.data[cell.hover, group.by, drop = TRUE]), collapse = '<br />')
} else {
NULL
}
}
)
}
# Run the thang
runGadget(app = ui, server = server)
}
#' @rdname LinkedPlots
#'
#' @aliases LinkedFeaturePlot
#'
#' @importFrom ggplot2 scale_fill_gradientn theme scale_alpha guides
#' scale_color_gradientn guide_colorbar
#'
#' @export
#' @concept visualization
#' @concept spatial
LinkedFeaturePlot <- function(
object,
feature,
dims = 1:2,
reduction = NULL,
image = NULL,
image.scale = "lowres",
slot = 'data',
alpha = c(0.1, 1),
combine = TRUE
) {
# Setup gadget UI
ui <- miniPage(
gadgetTitleBar(
title = 'LinkedFeaturePlot',
left = NULL
),
miniContentPanel(
fillRow(
plotOutput(
outputId = 'spatialplot',
height = '100%',
hover = hoverOpts(id = 'sphover', delay = 10, nullOutside = TRUE)
),
plotOutput(
outputId = 'dimplot',
height = '100%',
hover = hoverOpts(id = 'dimhover', delay = 10, nullOutside = TRUE)
),
height = '97%'
),
verbatimTextOutput(outputId = 'info')
)
)
# Prepare plotting data
cols <- SpatialColors(n = 100)
image <- image %||% DefaultImage(object = object)
cells.use <- Cells(x = object[[image]])
reduction <- reduction %||% DefaultDimReduc(object = object)
dims <- dims[1:2]
dims <- paste0(Key(object = object[[reduction]]), dims)
group.data <- FetchData(
object = object,
vars = feature,
cells = cells.use
)
coords <- GetTissueCoordinates(object = object[[image]], scale = image.scale)
embeddings <- Embeddings(object = object[[reduction]])[cells.use, dims]
sp_x <- colnames(coords)[1]
sp_y <- colnames(coords)[2]
dp_x <- dims[1]
dp_y <- dims[2]
# coordinates should be in image space, so need to flip y when setting or displaying info for points
flip_y <- function(pt) { if (!is.null(pt)) { pt$y <- max(coords[[sp_y]]) - (pt$y - min(coords[[sp_y]])) }; pt }
plot.data <- cbind(coords, group.data, embeddings)
# Setup the server
server <- function(input, output, session) {
plot.env <- reactiveValues()
# Handle events
observeEvent(
eventExpr = input$done,
handlerExpr = {
plots <- list(plot.env$spatialplot, plot.env$dimplot)
if (combine) {
plots <- wrap_plots(plots, ncol = 2)
}
stopApp(returnValue = plots)
}
)
# Set plots
output$spatialplot <- renderPlot(
expr = {
plot.env$spatialplot <- SingleSpatialPlot(
data = plot.data,
image = object[[image]],
col.by = feature,
pt.size.factor = 1.6,
crop = TRUE,
alpha.by = feature
) +
scale_fill_gradientn(name = feature, colours = cols) +
theme(legend.position = 'top') +
scale_alpha(range = alpha) +
guides(alpha = "none")
plot.env$spatialplot
}
)
output$dimplot <- renderPlot(
expr = {
plot.env$dimplot <- SingleDimPlot(
data = plot.data,
dims = dims,
col.by = feature
) +
scale_color_gradientn(name = feature, colours = cols, guide = 'colorbar') +
guides(color = guide_colorbar())
plot.env$dimplot
}
)
# Add hover text
output$info <- renderPrint(
expr = {
cell.hover <- rownames(x = nearPoints(
df = plot.data,
coordinfo = if (is.null(x = input[['sphover']])) {
input$dimhover
} else {
flip_y(input$sphover)
},
threshold = 10,
maxpoints = 1,
# specify plot-specific axis columns for nearPoints
xvar = if (is.null(x = input$sphover)) dp_x else sp_x,
yvar = if (is.null(x = input$sphover)) dp_y else sp_y
))
# TODO: Get newlines, extra information, and background color working
if (length(x = cell.hover) == 1) {
paste(cell.hover, paste('Expression:', plot.data[cell.hover, feature, drop = TRUE]), collapse = '<br />')
} else {
NULL
}
}
)
}
runGadget(app = ui, server = server)
}
#' Visualize clusters spatially and interactively
#'
#' @inheritParams SpatialPlot
#' @inheritParams DimPlot
#' @inheritParams LinkedPlots
#'
#' @return Returns final plot as a ggplot object
#'
#' @importFrom ggplot2 scale_alpha_ordinal
#' @importFrom miniUI miniPage miniButtonBlock miniTitleBarButton miniContentPanel
#' @importFrom shiny fillRow plotOutput verbatimTextOutput reactiveValues
#' observeEvent stopApp nearPoints renderPlot runGadget
#'
#' @export
#' @concept visualization
#' @concept spatial
#'
ISpatialDimPlot <- function(
object,
image = NULL,
image.scale = "lowres",
group.by = NULL,
alpha = c(0.3, 1)
) {
# Setup gadget UI
ui <- miniPage(
miniButtonBlock(miniTitleBarButton(
inputId = 'done',
label = 'Done',
primary = TRUE
)),
miniContentPanel(
fillRow(
plotOutput(
outputId = 'plot',
height = '100%',
click = clickOpts(id = 'click', clip = TRUE),
hover = hoverOpts(id = 'hover', delay = 10, nullOutside = TRUE)
),
height = '97%'
),
verbatimTextOutput(outputId = 'info')
)
)
# Get plotting data
# Prepare plotting data
image <- image %||% DefaultImage(object = object)
cells.use <- Cells(x = object[[image]])
group.by <- group.by %||% 'ident'
group.data <- FetchData(
object = object,
vars = group.by,
cells = cells.use
)
coords <- GetTissueCoordinates(object = object[[image]], scale = image.scale)
sp_x <- colnames(coords)[1]
sp_y <- colnames(coords)[2]
# coordinates should be in image space, so need to flip y when setting or displaying info for points
flip_y <- function(pt) { if (!is.null(pt)) { pt$y <- max(coords[[sp_y]]) - (pt$y - min(coords[[sp_y]])) }; pt }
scale.factor <- ScaleFactors(object[[image]])[[image.scale]]
plot.data <- cbind(coords, group.data)
plot.data$selected_ <- FALSE
Idents(object = object) <- group.by
# Set up the server
server <- function(input, output, session) {
click <- reactiveValues(pt = NULL)
plot.env <- reactiveValues(data = plot.data, alpha.by = NULL)
# Handle events
observeEvent(
eventExpr = input$done,
handlerExpr = stopApp(returnValue = plot.env$plot)
)
observeEvent(
eventExpr = input$click,
handlerExpr = {
click$pt <- flip_y(input$click)
clicked <- nearPoints(
df = plot.data,
coordinfo = click$pt,
threshold = 10,
maxpoints = 1,
xvar = sp_x,
yvar = sp_y
)
plot.env$data <- if (nrow(x = clicked) == 1) {
cell.clicked <- rownames(x = clicked)
group.clicked <- plot.data[cell.clicked, group.by, drop = TRUE]
idx.group <- which(x = plot.data[[group.by]] == group.clicked)
plot.data[idx.group, 'selected_'] <- TRUE
plot.data
} else {
plot.data
}
plot.env$alpha.by <- if (any(plot.env$data$selected_)) {
'selected_'
} else {
NULL
}
}
)
# Set plot
output$plot <- renderPlot(
expr = {
plot.env$plot <- SingleSpatialPlot(
data = plot.env$data,
image = object[[image]],
col.by = group.by,
crop = TRUE,
alpha.by = plot.env$alpha.by,
pt.size.factor = 1.6
) + scale_alpha_ordinal(range = alpha) + NoLegend()
plot.env$plot
}
)
# Add hover text
output$info <- renderPrint(
expr = {
hovered <- nearPoints(
df = plot.data,
coordinfo = flip_y(input$hover),
threshold = 10,
maxpoints = 1,
xvar = sp_x,
yvar = sp_y
)
if (nrow(hovered) == 1) {
cell.hover <- rownames(hovered)
# SingleSpatialPlot relies on the spatial coordinates appearing
# in the first two columns of the returned data.frame - it's kinda
# fragile but we're obligated to use the same behaviour here
coords.hover <- hovered[1, colnames(coords)[1:2]] / scale.factor
group.hover <- hovered[1, group.by]
sprintf(
"Cell: %s, Group: %s, Coordinates: (%.2f, %.2f)",
cell.hover,
group.hover,
coords.hover[[1]],
coords.hover[[2]]
)
} else {
NULL
}
}
)
}
runGadget(app = ui, server = server)
}
#' Visualize features spatially and interactively
#'
#' @inheritParams SpatialPlot
#' @inheritParams FeaturePlot
#' @inheritParams LinkedPlots
#'
#' @return Returns final plot as a ggplot object
#'
#' @importFrom ggplot2 scale_fill_gradientn theme scale_alpha guides
#' @importFrom miniUI miniPage miniButtonBlock miniTitleBarButton miniContentPanel
#' @importFrom shiny fillRow sidebarPanel sliderInput selectInput reactiveValues
#' observeEvent stopApp observe updateSelectInput plotOutput renderPlot runGadget
#'
#' @export
#' @concept visualization
#' @concept spatial
ISpatialFeaturePlot <- function(
object,
feature,
image = NULL,
image.scale = "lowres",
slot = 'data',
alpha = c(0.1, 1)
) {
# Set inital data values
assay.keys <- Key(object = object)[Assays(object = object)]
keyed <- sapply(X = assay.keys, FUN = grepl, x = feature)
assay <- if (any(keyed)) {
names(x = which(x = keyed))[1]
} else {
DefaultAssay(object = object)
}
features <- sort(x = rownames(x = GetAssayData(
object = object,
layer = slot,
assay = assay
)))
feature.label <- 'Feature to visualize'
assays.use <- vapply(
X = Assays(object = object),
FUN = function(x) {
return(!IsMatrixEmpty(x = GetAssayData(
object = object,
layer = slot,
assay = x
)))
},
FUN.VALUE = logical(length = 1L)
)
assays.use <- sort(x = Assays(object = object)[assays.use])
# Setup gadget UI
ui <- miniPage(
miniButtonBlock(miniTitleBarButton(
inputId = 'done',
label = 'Done',
primary = TRUE
)),
miniContentPanel(
fillRow(
sidebarPanel(
sliderInput(
inputId = 'alpha',
label = 'Alpha intensity',
min = 0,
max = max(alpha),
value = min(alpha),
step = 0.01,
width = '100%'
),
sliderInput(
inputId = 'pt.size',
label = 'Point size',
min = 0,
max = 5,
value = 1.6,
step = 0.1,
width = '100%'
),
selectInput(
inputId = 'assay',
label = 'Assay',
choices = assays.use,
selected = assay,
selectize = FALSE,
width = '100%'
),
selectInput(
inputId = 'feature',
label = feature.label,
choices = features,
selected = feature,
selectize = FALSE,
width = '100%'
),
selectInput(
inputId = 'palette',
label = 'Color scheme',
choices = names(x = FeaturePalettes),
selected = 'Spatial',
selectize = FALSE,
width = '100%'
),
width = '100%'
),
plotOutput(outputId = 'plot', height = '100%'),
flex = c(1, 4)
)
)
)
# Prepare plotting data
image <- image %||% DefaultImage(object = object)
cells.use <- Cells(x = object[[image]])
coords <- GetTissueCoordinates(object = object[[image]], scale = image.scale)
feature.data <- FetchData(
object = object,
vars = feature,
cells = cells.use,
layer = slot
)
plot.data <- cbind(coords, feature.data)
server <- function(input, output, session) {
plot.env <- reactiveValues(
data = plot.data,
feature = feature,
palette = 'Spatial'
)
# Observe events
observeEvent(
eventExpr = input$done,
handlerExpr = stopApp(returnValue = plot.env$plot)
)
observe(x = {
assay <- input$assay
feature.use <- input$feature
features.assay <- sort(x = rownames(x = GetAssayData(
object = object,
layer = slot,
assay = assay
)))
feature.use <- ifelse(
test = feature.use %in% features.assay,
yes = feature.use,
no = features.assay[1]
)
updateSelectInput(
session = session,
inputId = 'assay',
label = 'Assay',
choices = assays.use,
selected = assay
)
updateSelectInput(
session = session,
inputId = 'feature',
label = feature.label,
choices = features.assay,
selected = feature.use
)
})
observe(x = {
feature.use <- input$feature
try(
expr = {
feature.data <- FetchData(
object = object,
vars = paste0(Key(object = object[[input$assay]]), feature.use),
cells = cells.use,
layer = slot
)
colnames(x = feature.data) <- feature.use
plot.env$data <- cbind(coords, feature.data)
plot.env$feature <- feature.use
},
silent = TRUE
)
})
observe(x = {
plot.env$palette <- input$palette
})
# Create plot
output$plot <- renderPlot(expr = {
plot.env$plot <- SingleSpatialPlot(
data = plot.env$data,
image = object[[image]],
col.by = plot.env$feature,
pt.size.factor = input$pt.size,
crop = TRUE,
alpha.by = plot.env$feature
) +
# scale_fill_gradientn(name = plot.env$feature, colours = cols) +
scale_fill_gradientn(name = plot.env$feature, colours = FeaturePalettes[[plot.env$palette]]) +
theme(legend.position = 'top') +
scale_alpha(range = c(input$alpha, 1)) +
guides(alpha = "none")
plot.env$plot
})
}
runGadget(app = ui, server = server)
}
#' Interactive Spatial Cell Selection Tool
#'
#' Launch an interactive gadget for lasso-based cell selection from a spatial Seurat object.
#' Supports Visium, SlideSeq, and Vizgen data. Returns the cell names of the selected subset,
#' suitable for downstream subsetting or analysis.
#'
#' @note This function requires the
#' \href{https://cran.r-project.org/package=plotly}{\pkg{plotly}},
#' \href{https://cran.r-project.org/package=magrittr}{\pkg{magrittr}},
#' and \href{https://cran.r-project.org/package=base64enc}{\pkg{base64enc}} packages
#' to be installed. It also requires \pkg{shiny} and \pkg{miniUI} for the interactive UI.
#'
#' @param object A \code{\link[SeuratObject]{Seurat}} object with spatial data.
#' @param image Name of the spatial image stored in the object. If \code{NULL}, uses the default image for the object.
#' @param image.scale Character. Which image scaling factor to use for spatial coordinate transformation (\code{"lowres"} by default).
#' @param group.by Metadata variable (column name) to use for coloring cell points (e.g., cluster assignment). If \code{NULL}, uses \code{"seurat_clusters"} if available, otherwise all cells are grouped together.
#' @param alpha Numeric transparency value for cell points (default \code{1.0}).
#' @param pt.size.factor Numeric scaling factor for point size (default \code{1.0}).
#' @param overlay_image Logical; if \code{TRUE}, overlays the tissue image in the background of the plot (default \code{TRUE}).
#'
#' @importFrom grDevices png dev.off
#' @importFrom miniUI miniPage gadgetTitleBar miniContentPanel
#' @importFrom shiny uiOutput reactiveVal renderUI tags observeEvent stopApp runGadget
#'
#' @return A character vector of cell names selected via lasso, which can be used to subset the object.
#' @export
#'
#' @examples
#' \dontrun{
#' selected_cells <- InteractiveSpatialPlot(object = brain)
#' selected_cells <- InteractiveSpatialPlot(object = brain, overlay_image = FALSE)
#' }
InteractiveSpatialPlot <- function(
object,
image = NULL,
image.scale = "lowres",
group.by = NULL,
alpha = 1.0,
pt.size.factor = 1.0,
overlay_image = TRUE
) {
# Check for required packages, stop with clear message if missing
required_pkgs <- c("plotly", "magrittr", "base64enc", "shiny")
missing_pkgs <- required_pkgs[
!vapply(required_pkgs, requireNamespace, quietly = TRUE, FUN.VALUE = logical(1))
]
if (length(missing_pkgs) > 0) {
stop(
"InteractiveSpatialPlot() functionality requires these packages to be installed: ",
paste0("'", missing_pkgs, "'", collapse = ", "),
call. = FALSE
)
}
# Import magrittr pipe locally
`%>%` <- magrittr::`%>%`
# Use provided image name or fallback to default
image <- image %||% DefaultImage(object)
# Sanity check: requested image must exist in the object
if (!image %in% names(object@images)) {
stop("Image '", image, "' not found. Available image(s): ", paste(names(object@images), collapse = ", "))
}
# Retrieve the spatial image object
image_obj <- object[[image]]
# Determine image technology type (Visium, SlideSeq, or Vizgen)
img_class <- class(image_obj)[1]
if (img_class %in% c("VisiumV1", "VisiumV2")) {
type <- "visium"
} else if (img_class == "SlideSeq") {
type <- "slideseq"
} else if (img_class == "FOV") {
type <- "vizgen"
} else {
stop("Unrecognized image class: ", img_class)
}
# Extract and scale cell coordinates according to image type
if (type == "visium") {
# For Visium: coordinates stored in centroids, need scaling
if (!"boundaries" %in% slotNames(image_obj)) {
stop("Image object does not have a 'boundaries' slot; check if data is truly Visium data")
}
centroids <- image_obj@boundaries$centroids
coords <- setNames(as.data.frame(centroids@coords), c("x", "y"))
coords$cell <- centroids@cells
# Scale coordinates to match image pixel units
scale.factor <- Seurat::ScaleFactors(image_obj)[[image.scale]]
if (is.null(scale.factor)) stop("Scale factor for '", image.scale, "' not found")
coords$x_raw <- coords$x # Store original, unscaled x
coords$y_raw <- coords$y # Store original, unscaled y
coords$x <- coords$x * scale.factor
coords$y <- coords$y * scale.factor
} else if (type == "slideseq") {
# For Slide-seq: coordinates are stored directly
if (!"coordinates" %in% slotNames(image_obj)) {
stop("Image object does not have a 'coordinates' slot; check if data is truly Slide-seq data")
}
coords <- as.data.frame(image_obj@coordinates)
coords$cell <- rownames(coords)
colnames(coords)[1:2] <- c("x", "y")
coords$x_raw <- coords$x
coords$y_raw <- coords$y
} else if (type == "vizgen") {
# For Vizgen: coordinates in centroids
if (!"boundaries" %in% slotNames(image_obj)) {
stop("Vizgen FOV missing 'boundaries' slot")
}
centroids <- image_obj@boundaries$centroids
coords <- as.data.frame(centroids@coords)
colnames(coords) <- c("x", "y")
coords$cell <- centroids@cells
coords$x_raw <- coords$x
coords$y_raw <- coords$y
}
# Get cell-level metadata for grouping/labeling
meta <- [email hidden]
# If group.by not given, use 'seurat_clusters' if available; otherwise group all together
if (is.null(group.by)) {
group.by <- if ("seurat_clusters" %in% colnames(meta)) "seurat_clusters" else "all"
}
# Assign group/cluster for coloring the plot
if (group.by != "all") {
coords$group <- meta[coords$cell, group.by]
} else {
coords$group <- "all"
}
# Compose hover text: show cell name, original (x, y) coordinates, rounded for clarity
coords$hover <- paste0(
"Cell: ", coords$cell,
"<br>x: ", round(coords$x_raw, 1),
", y: ", round(coords$y_raw, 1)
)
# Prepare background tissue image as a base64-encoded PNG (if available and enabled)
base64_image <- NULL
img_width <- NULL
img_height <- NULL
if (overlay_image) {
# Only attempt to overlay image if compatible type and slot present
if (type == "visium" && "image" %in% slotNames(image_obj)) {
img_raster <- image_obj@image
} else if (
type == "vizgen" &&
"boundaries" %in% slotNames(image_obj) &&
"centroids" %in% slotNames(image_obj@boundaries) &&
"image" %in% slotNames(image_obj@boundaries$centroids)
) {
img_raster <- image_obj@boundaries$centroids@image
}
# Convert the raster image array to base64 PNG (for embedding in plotly)
if (exists("img_raster")) {
img_width <- dim(img_raster)[2]
img_height <- dim(img_raster)[1]
temp_png <- tempfile(fileext = ".png")
png(temp_png, width = img_width, height = img_height)
grid::grid.raster(img_raster)
dev.off()
img_bytes <- readBin(temp_png, "raw", file.info(temp_png)$size)
base64_image <- paste0("data:image/png;base64,", base64enc::base64encode(img_bytes))
}
}
# Calculate custom axis tick positions and labels to show original coordinates
# This is necessary as points are downscaled to fit on the tissue image
# However, to best retain their original spatial orientation, we plot
# the original coordinate scale on the axis
create_axis_ticks <- function(scaled_coords, raw_coords, n_ticks = 6) {
# Get range of scaled and raw coordinates
scaled_range <- range(scaled_coords, na.rm = TRUE)
raw_range <- range(raw_coords, na.rm = TRUE)
# Create tick positions in the raw coordinate space
raw_ticks <- pretty(raw_range, n = n_ticks)
# Calculate corresponding scaled positions
# Linear interpolation from raw to scaled coordinates
scale_factor <- diff(scaled_range) / diff(raw_range)
scaled_ticks <- (raw_ticks - raw_range[1]) * scale_factor + scaled_range[1]
return(list(tickvals = scaled_ticks, ticktext = as.character(raw_ticks)))
}
# Create custom axis ticks for both x and y axes
x_ticks <- create_axis_ticks(coords$x, coords$x_raw)
y_ticks <- create_axis_ticks(coords$y, coords$y_raw)
# Set up the gadget UI with a plotly output area
ui <- miniPage(
gadgetTitleBar("Select a subset of cells"),
miniContentPanel(
plotly::plotlyOutput("plot", height = "100%"),
shiny::tags$div(
shiny::uiOutput("selection_count"),
style = "position:absolute; bottom:8px; right:10px; padding:4px 6px; background:rgba(255,255,255,0.8); font-size:12px; border-radius:3px; pointer-events:none;"
)
)
)
# Shiny gadget server logic for interactive plot and lasso selection
server <- function(input, output, session) {
current_selection <- shiny::reactiveVal(coords$cell)
# Render the interactive plotly scattergl plot
output$plot <- plotly::renderPlotly({
plt <- plotly::plot_ly(
data = coords,
x = ~x,
y = ~y,
color = ~group, # Color by group/cluster if available
key = ~cell, # Store cell names for selection retrieval
type = "scattergl", # Use WebGL for performance with large datasets
mode = "markers",
marker = list(size = 2 * pt.size.factor), # Default pt size is 2
text = ~hover, # Show hover info (cellid + coordinates)
hoverinfo = "text",
alpha = alpha # Global transparency
)
# Overlay the tissue image as background if available
if (!is.null(base64_image)) {
plt <- plt %>% plotly::layout(
images = list(
list(
source = base64_image,
xref = "x", yref = "y", # Anchor to data coordinates
x = 0,
y = 0,
sizex = img_width,
sizey = img_height,
sizing = "stretch",
opacity = 0.6,
layer = "below"
)
)
)
}
# Lock axes to same scale and reverse y for image alignment
# Set lasso mode and custom axis labels
plt <- plt %>% plotly::layout(
dragmode = "lasso",
yaxis = list(
autorange = "reversed",
scaleanchor = "x",
title = "y",
tickvals = y_ticks$tickvals,
ticktext = y_ticks$ticktext
),
xaxis = list(
scaleanchor = "y",
title = "x",
tickvals = x_ticks$tickvals,
ticktext = x_ticks$ticktext
)
)
plt
})
observeEvent(plotly::event_data("plotly_selected"), {
selected <- plotly::event_data("plotly_selected")
if (is.null(selected) || NROW(selected) == 0) {
current_selection(NULL)
} else {
keys <- selected$key
keys <- keys[!is.na(keys)]
current_selection(keys)
}
}, ignoreInit = TRUE)
output$selection_count <- shiny::renderUI({
shiny::tags$span(paste0("Selected cells: ", NROW(current_selection())))
})
# When user clicks "Done", retrieve lasso selection and close gadget
observeEvent(input$done, {
stopApp(current_selection())
})
# When user clicks "Cancel", exit gadget and return NULL
observeEvent(input$cancel, {
stopApp(NULL)
})
}
# Launch the interactive gadget
runGadget(ui, server)
}
#' Visualize spatial clustering and expression data.
#'
#' SpatialPlot plots a feature or discrete grouping (e.g. cluster assignments) as
#' spots over the image that was collected. We also provide SpatialFeaturePlot
#' and SpatialDimPlot as wrapper functions around SpatialPlot for a consistent
#' naming framework.
#'
#' @inheritParams HoverLocator
#' @param object A Seurat object
#' @param group.by Name of meta.data column to group the data by
#' @param features Name of the feature to visualize. Provide either group.by OR
#' features, not both.
#' @param images Name of the images to use in the plot(s)
#' @param cols Vector of colors, each color corresponds to an identity class.
#' This may also be a single character or numeric value corresponding to a
#' palette as specified by \code{\link[RColorBrewer]{brewer.pal.info}}. By
#' default, ggplot2 assigns colors
#' @param image.alpha Adjust the opacity of the background images. Set to 0 to
#' remove.
#' @param image.scale Choose the scale factor ("lowres"/"hires") to apply in
#' order to matchthe plot with the specified `image` - defaults to "lowres"
#' @param crop Crop the plot in to focus on points plotted. Set to \code{FALSE} to show
#' entire background image.
#' @param slot If plotting a feature, which data slot to pull from (counts,
#' data, or scale.data)
#' @param keep.scale How to handle the color scale across multiple plots. Options are:
#' \itemize{
#' \item \dQuote{feature} (default; by row/feature scaling): The plots for
#' each individual feature are scaled to the maximum expression of the
#' feature across the conditions provided to \code{split.by}
#' \item \dQuote{all} (universal scaling): The plots for all features and
#' conditions are scaled to the maximum expression value for the feature
#' with the highest overall expression
#' \item \code{NULL} (no scaling): Each individual plot is scaled to the
#' maximum expression value of the feature in the condition provided to
#' \code{split.by}; be aware setting \code{NULL} will result in color
#' scales that are not comparable between plots
#' }
#' @param min.cutoff,max.cutoff Vector of minimum and maximum cutoff
#' values for each feature, may specify quantile in the form of 'q##' where '##'
#' is the quantile (eg, 'q1', 'q10')
#' @param cells.highlight A list of character or numeric vectors of cells to
#' highlight. If only one group of cells desired, can simply pass a vector
#' instead of a list. If set, colors selected cells to the color(s) in
#' cols.highlight
#' @param cols.highlight A vector of colors to highlight the cells as; ordered
#' the same as the groups in cells.highlight; last color corresponds to
#' unselected cells.
#' @param facet.highlight When highlighting certain groups of cells, split each
#' group into its own plot
#' @param label Whether to label the clusters
#' @param label.size Sets the size of the labels
#' @param label.color Sets the color of the label text
#' @param label.box Whether to put a box around the label text (geom_text vs
#' geom_label)
#' @param repel Repels the labels to prevent overlap
#' @param ncol Number of columns if plotting multiple plots
#' @param combine Combine plots into a single gg object; note that if TRUE;
#' themeing will not work when plotting multiple features/groupings
#' @param pt.size.factor Scale the size of the spots.
#' @param alpha Controls opacity of spots. Provide as a vector specifying the
#' min and max for SpatialFeaturePlot. For SpatialDimPlot, provide a single
#' alpha value for each plot.
#' @param shape Control the shape of the spots - same as the ggplot2 parameter.
#' The default is 21, which plots circles - use 22 to plot squares.
#' @param stroke Control the width of the border around the spots
#' @param stroke.alpha Control the opacity of spot borders (when stroke is specified).
#' Set to \code{NA} to use the same alpha as the fill.
#' @param interactive Launch an interactive SpatialDimPlot or SpatialFeaturePlot
#' session, see \code{\link{ISpatialDimPlot}} or
#' \code{\link{ISpatialFeaturePlot}} for more details
#' @param do.identify,do.hover DEPRECATED in favor of \code{interactive}
#' @param identify.ident DEPRECATED
#' @param plot_segmentations Define whether plot should plot centroids or segmentations
#'
#' @return If \code{do.identify}, either a vector of cells selected or the object
#' with selected cells set to the value of \code{identify.ident} (if set). Else,
#' if \code{do.hover}, a plotly object with interactive graphics. Else, a ggplot
#' object
#'
#' @importFrom ggplot2 scale_fill_gradientn ggtitle theme element_text scale_alpha
#' @importFrom patchwork wrap_plots
#' @export
#' @concept visualization
#' @concept spatial
#'
#' @examples
#' \dontrun{
#' # For functionality analagous to FeaturePlot
#' SpatialPlot(seurat.object, features = "MS4A1")
#' SpatialFeaturePlot(seurat.object, features = "MS4A1")
#'
#' # For functionality analagous to DimPlot
#' SpatialPlot(seurat.object, group.by = "clusters")
#' SpatialDimPlot(seurat.object, group.by = "clusters")
#' }
#'
SpatialPlot <- function(
object,
group.by = NULL,
features = NULL,
images = NULL,
cols = NULL,
image.alpha = 1,
image.scale = "lowres",
crop = TRUE,
slot = 'data',
keep.scale = "feature",
min.cutoff = NA,
max.cutoff = NA,
cells.highlight = NULL,
cols.highlight = c('#DE2D26', 'grey50'),
facet.highlight = FALSE,
label = FALSE,
label.size = 5,
label.color = 'white',
label.box = TRUE,
repel = FALSE,
ncol = NULL,
combine = TRUE,
pt.size.factor = 1.6,
alpha = c(1, 1),
shape = 21,
stroke = NA,
stroke.alpha = NA,
interactive = FALSE,
do.identify = FALSE,
identify.ident = NULL,
do.hover = FALSE,
information = NULL,
plot_segmentations = FALSE
) {
if (isTRUE(x = do.hover) || isTRUE(x = do.identify)) {
warning(
"'do.hover' and 'do.identify' are deprecated as we are removing plotly-based interactive graphics, use 'interactive' instead for Shiny-based interactivity",
call. = FALSE,
immediate. = TRUE
)
interactive <- TRUE
}
if (!is.null(x = group.by) & !is.null(x = features)) {
stop("Please specific either group.by or features, not both.")
}
images <- images %||% Images(object = object, assay = DefaultAssay(object = object))
if (length(x = images) == 0) {
images <- Images(object = object)
}
if (length(x = images) < 1) {
stop("Could not find any spatial image information")
}
# Check keep.scale param for valid entries
if (!(is.null(x = keep.scale)) && !(keep.scale %in% c("feature", "all"))) {
stop("`keep.scale` must be set to either `feature`, `all`, or NULL")
}
cells <- unique(CellsByImage(object, images = images, unlist = TRUE))
if (is.null(x = features)) {
if (interactive) {
# default alpha is 1 but interactive plotting requires
# a range for proper cluster selection highlighting
if (identical(alpha, c(1, 1))) {
alpha <- c(0.1, 1)
}
tryCatch(
expr = {
return(ISpatialDimPlot(
object = object,
image = images[1],
image.scale = image.scale,
group.by = group.by,
alpha = alpha
))
},
error = function(e) {
# error can occur when image and assay don't match
# or when the default assay set doesn't have data corresponding to the default ident etc.
if (grepl("arguments imply differing number of rows", conditionMessage(e))) {
stop(
"Cells were removed due to missing data; check if the specified image and assay are correct.\n",
call. = FALSE
)
} else {
stop(e)
}
}
)
}
group.by <- group.by %||% 'ident'
object[['ident']] <- Idents(object = object)
data <- object[[group.by]]
data <- data[cells,,drop=F]
for (group in group.by) {
if (!is.factor(x = data[, group])) {
data[, group] <- factor(x = data[, group])
}
}
} else {
if (interactive) {
return(ISpatialFeaturePlot(
object = object,
feature = features[1],
image = images[1],
image.scale = image.scale,
slot = slot,
alpha = alpha
))
}
data <- FetchData(
object = object,
vars = features,
cells = cells,
layer = slot,
clean = FALSE
)
features <- colnames(x = data)
# Determine cutoffs
min.cutoff <- mapply(
FUN = function(cutoff, feature) {
return(ifelse(
test = is.na(x = cutoff),
yes = min(data[, feature]),
no = cutoff
))
},
cutoff = min.cutoff,
feature = features
)
max.cutoff <- mapply(
FUN = function(cutoff, feature) {
return(ifelse(
test = is.na(x = cutoff),
yes = max(data[, feature]),
no = cutoff
))
},
cutoff = max.cutoff,
feature = features
)
check.lengths <- unique(x = vapply(
X = list(features, min.cutoff, max.cutoff),
FUN = length,
FUN.VALUE = numeric(length = 1)
))
if (length(x = check.lengths) != 1) {
stop("There must be the same number of minimum and maximum cuttoffs as there are features")
}
# Apply cutoffs
data <- sapply(
X = 1:ncol(x = data),
FUN = function(index) {
data.feature <- as.vector(x = data[, index])
min.use <- SetQuantile(cutoff = min.cutoff[index], data.feature)
max.use <- SetQuantile(cutoff = max.cutoff[index], data.feature)
data.feature[data.feature < min.use] <- min.use
data.feature[data.feature > max.use] <- max.use
return(data.feature)
}
)
colnames(x = data) <- features
rownames(x = data) <- cells
}
features <- colnames(x = data)
colnames(x = data) <- features
rownames(x = data) <- cells
facet.highlight <- facet.highlight && (!is.null(x = cells.highlight) && is.list(x = cells.highlight))
if (do.hover) {
if (length(x = images) > 1) {
images <- images[1]
warning(
"'do.hover' requires only one image, using image ",
images,
call. = FALSE,
immediate. = TRUE
)
}
if (length(x = features) > 1) {
features <- features[1]
type <- ifelse(test = is.null(x = group.by), yes = 'feature', no = 'grouping')
warning(
"'do.hover' requires only one ",
type,
", using ",
features,
call. = FALSE,
immediate. = TRUE
)
}
if (facet.highlight) {
warning(
"'do.hover' requires no faceting highlighted cells",
call. = FALSE,
immediate. = TRUE
)
facet.highlight <- FALSE
}
}
if (facet.highlight) {
if (length(x = images) > 1) {
images <- images[1]
warning(
"Faceting the highlight only works with a single image, using image ",
images,
call. = FALSE,
immediate. = TRUE
)
}
ncols <- length(x = cells.highlight)
} else {
ncols <- length(x = images)
}
plots <- vector(
mode = "list",
length = length(x = features) * ncols
)
# Get max across all features
if (!(is.null(x = keep.scale)) && keep.scale == "all") {
max.feature.value <- max(apply(data, 2, function(x) max(x, na.rm = TRUE)))
}
for (i in 1:ncols) {
plot.idx <- i
image.idx <- ifelse(test = facet.highlight, yes = 1, no = i)
image.use <- object[[images[[image.idx]]]]
is_visium_v2 <- inherits(image.use, "VisiumV2")
old_axis_orientation <- (!.hasSlot(image.use, "coords_x_orientation")) || (.hasSlot(image.use, "coords_x_orientation") && (slot(image.use, "coords_x_orientation") != 'horizontal'))
if (is_visium_v2 && old_axis_orientation) {
stop(
"Please run `UpdateSeuratObject` on your Seurat object first to ensure that data aligns to the image ", images[[image.idx]], " when plotting.",
call. = TRUE
)
}
# When plotting segmentations, set the default boundary (temporarily) to segmentations
if (plot_segmentations == TRUE && inherits(image.use, "VisiumV2") &&
"segmentations" %in% names(image.use)) {
db <- DefaultBoundary(image.use)
on.exit(DefaultBoundary(image.use) <- db, add = TRUE) # Reset on exit
DefaultBoundary(image.use) <- "segmentations"
}
coordinates <- GetTissueCoordinates(object = image.use, scale = image.scale)
highlight.use <- if (facet.highlight) {
cells.highlight[i]
} else {
cells.highlight
}
for (j in seq_along(features)) {
cols.unset <- is.factor(x = data[, features[j]]) && is.null(x = cols)
if (cols.unset) {
cols <- hue_pal()(n = length(x = levels(x = data[, features[j]])))
names(x = cols) <- levels(x = data[, features[j]])
}
# Get feature max for individual feature
if (!(is.null(x = keep.scale)) && keep.scale == "feature" && !inherits(x = data[, features[j]], what = "factor") ) {
max.feature.value <- max(data[, features[j]])
}
# Check if object is of type Visium and contains segmentations
has_visium_segm_data <- inherits(image.use, "VisiumV2") &&
!is.null(image.use@boundaries$segmentations) &&
"sf.data" %in% slotNames(image.use@boundaries$segmentations)
# GetTissueCoordinates will not always return a "cell" column (e.g., Visium V1)
if (!("cell" %in% colnames(x = coordinates))) {
coordinates$cell <- rownames(x = coordinates)
}
idx <- match(coordinates$cell, rownames(x = data))
plot.data <- cbind(coordinates, data[idx, features[j], drop = FALSE])
plot <- SingleSpatialPlot(
data = plot.data,
image = image.use,
image.scale = image.scale,
image.alpha = image.alpha,
col.by = features[j],
cols = cols,
alpha.by = if (is.null(x = group.by)) {
features[j]
} else {
NULL
},
pt.alpha = if (!is.null(x = group.by)) {
alpha[j]
} else {
NULL
},
geom = if (inherits(x = image.use, what = "STARmap")) {
"poly_starmap"
} else if (has_visium_segm_data && plot_segmentations) {
"poly"
} else {
"spatial"
},
cells.highlight = highlight.use,
cols.highlight = cols.highlight,
pt.size.factor = pt.size.factor,
shape = shape,
stroke = stroke,
stroke.alpha = stroke.alpha,
crop = crop
)
if (is.null(x = group.by)) {
plot <- plot +
scale_fill_gradientn(
name = features[j],
colours = SpatialColors(n = 100)
) +
theme(legend.position = 'top') +
scale_alpha(range = alpha) +
guides(alpha = "none")
} else if (label) {
plot <- LabelClusters(
plot = plot,
id = ifelse(
test = is.null(x = cells.highlight),
yes = features[j],
no = 'highlight'
),
geom = if (inherits(x = image.use, what = "STARmap") || (has_visium_segm_data && plot_segmentations)) {
'GeomPolygon'
} else {
'GeomSpatial'
},
repel = repel,
size = label.size,
color = label.color,
box = label.box,
position = "nearest"
)
}
if (j == 1 && length(x = images) > 1 && !facet.highlight) {
plot <- plot +
ggtitle(label = images[[image.idx]]) +
theme(plot.title = element_text(hjust = 0.5))
}
if (facet.highlight) {
plot <- plot +
ggtitle(label = names(x = cells.highlight)[i]) +
theme(plot.title = element_text(hjust = 0.5)) +
NoLegend()
}
if (has_visium_segm_data && plot_segmentations && !is.null(group.by)) {
# Add legend guides to show filled squares next to labels when plotting segmentations
plot <- plot + guides(fill = guide_legend(override.aes = list(alpha = 1, color = "black", linewidth = 0.2, size = 2)))
}
# Plot multiple images depending on keep.scale
if (!(is.null(x = keep.scale)) && !inherits(x = data[, features[j]], "factor")) {
plot <- suppressMessages(plot & scale_fill_gradientn(colors = SpatialColors(n = 100), limits = c(NA, max.feature.value)))
}
plots[[plot.idx]] <- plot
plot.idx <- plot.idx + ncols
if (cols.unset) {
cols <- NULL
}
}
}
if (combine) {
if (!is.null(x = ncol)) {
return(wrap_plots(plots = plots, ncol = ncol))
}
if (length(x = images) > 1) {
return(wrap_plots(plots = plots, ncol = length(x = images)))
}
return(wrap_plots(plots = plots))
}
return(plots)
}
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Other plotting functions
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#' Plot the Barcode Distribution and Calculated Inflection Points
#'
#' This function plots the calculated inflection points derived from the barcode-rank
#' distribution.
#'
#' See [CalculateBarcodeInflections()] to calculate inflection points and
#' [SubsetByBarcodeInflections()] to subsequently subset the Seurat object.
#'
#' @param object Seurat object
#'
#' @return Returns a `ggplot2` object showing the by-group inflection points and provided
#' (or default) rank threshold values in grey.
#'
#' @importFrom methods slot
#' @importFrom cowplot theme_cowplot
#' @importFrom ggplot2 ggplot geom_line geom_vline
#'
#' @export
#' @concept visualization
#'
#' @author Robert A. Amezquita, \email{[email hidden]}
#' @seealso \code{\link{CalculateBarcodeInflections}} \code{\link{SubsetByBarcodeInflections}}
#'
#' @examples
#' data("pbmc_small")
#' pbmc_small <- CalculateBarcodeInflections(pbmc_small, group.column = 'groups')
#' BarcodeInflectionsPlot(pbmc_small)
#'
BarcodeInflectionsPlot <- function(object) {
cbi.data <- Tool(object = object, slot = 'CalculateBarcodeInflections')
if (is.null(x = cbi.data)) {
stop("Barcode inflections not calculated, please run CalculateBarcodeInflections")
}
## Extract necessary data frames
inflection_points <- cbi.data$inflection_points
barcode_distribution <- cbi.data$barcode_distribution
threshold_values <- cbi.data$threshold_values
# Set a cap to max rank to avoid plot being overextended
if (threshold_values$rank[[2]] > max(barcode_distribution$rank, na.rm = TRUE)) {
threshold_values$rank[[2]] <- max(barcode_distribution$rank, na.rm = TRUE)
}
## Infer the grouping/barcode variables
group_var <- colnames(x = barcode_distribution)[1]
barcode_var <- colnames(x = barcode_distribution)[2]
barcode_distribution[, barcode_var] <- log10(x = barcode_distribution[, barcode_var] + 1)
## Make the plot
plot <- ggplot(
data = barcode_distribution,
mapping = aes(
x = .data[['rank']],
y = .data[[barcode_var]],
group = .data[[group_var]],
colour = .data[[group_var]]
)
) +
geom_line() +
geom_vline(
data = threshold_values,
aes(xintercept = .data[['rank']]),
linetype = "dashed",
colour = 'grey60',
size = 0.5
) +
geom_vline(
data = inflection_points,
mapping = aes(
xintercept = .data[['rank']],
group = .data[[group_var]],
colour = .data[[group_var]]
),
linetype = "dashed"
) +
theme_cowplot()
return(plot)
}
#' Dot plot visualization
#'
#' Intuitive way of visualizing how feature expression changes across different
#' identity classes (clusters). The size of the dot encodes the percentage of
#' cells within a class, while the color encodes the AverageExpression level
#' across all cells within a class (blue is high).
#'
#' @param object Seurat object
#' @param assay Name of assay to use, defaults to the active assay
#' @param features Input vector of features, or named list of feature vectors
#' if feature-grouped panels are desired (replicates the functionality of the
#' old SplitDotPlotGG)
#' @param cols Colors to plot: the name of a palette from
#' \code{RColorBrewer::brewer.pal.info}, a pair of colors defining a gradient,
#' or 3+ colors defining multiple gradients (if split.by is set)
#' @param col.min Minimum scaled average expression threshold (everything
#' smaller will be set to this)
#' @param col.max Maximum scaled average expression threshold (everything larger
#' will be set to this)
#' @param dot.min The fraction of cells at which to draw the smallest dot
#' (default is 0). All cell groups with less than this expressing the given
#' gene will have no dot drawn.
#' @param dot.scale Scale the size of the points, similar to cex
#' @param idents Identity classes to include in plot (default is all)
#' @param group.by Factor to group the cells by
#' @param split.by A factor in object metadata to split the plot by, pass 'ident'
#' to split by cell identity
#' see \code{\link{FetchData}} for more details
#' @param cluster.idents Whether to order identities by hierarchical clusters
#' based on given features, default is FALSE
#' @param scale Determine whether the data is scaled, TRUE for default
#' @param scale.by Scale the size of the points by 'size' or by 'radius'
#' @param scale.min Set lower limit for scaling, use NA for default
#' @param scale.max Set upper limit for scaling, use NA for default
#'
#' @return A ggplot object
#'
#' @importFrom grDevices colorRampPalette
#' @importFrom cowplot theme_cowplot
#' @importFrom ggplot2 ggplot geom_point scale_size scale_radius
#' theme element_blank labs scale_color_identity scale_color_distiller
#' scale_color_gradient guides guide_legend guide_colorbar
#' facet_grid unit
#' @importFrom scattermore geom_scattermore
#' @importFrom stats dist hclust
#' @importFrom RColorBrewer brewer.pal.info
#'
#' @export
#' @concept visualization
#'
#' @aliases SplitDotPlotGG
#' @seealso \code{RColorBrewer::brewer.pal.info}
#'
#' @examples
#' data("pbmc_small")
#' cd_genes <- c("CD247", "CD3E", "CD9")
#' DotPlot(object = pbmc_small, features = cd_genes)
#' pbmc_small[['groups']] <- sample(x = c('g1', 'g2'), size = ncol(x = pbmc_small), replace = TRUE)
#' DotPlot(object = pbmc_small, features = cd_genes, split.by = 'groups')
#'
DotPlot <- function(
object,
features,
assay = NULL,
cols = c("lightgrey", "blue"),
col.min = -2.5,
col.max = 2.5,
dot.min = 0,
dot.scale = 6,
idents = NULL,
group.by = NULL,
split.by = NULL,
cluster.idents = FALSE,
scale = TRUE,
scale.by = 'radius',
scale.min = NA,
scale.max = NA
) {
assay <- assay %||% DefaultAssay(object = object)
DefaultAssay(object = object) <- assay
split.colors <- !is.null(x = split.by) && !any(cols %in% rownames(x = brewer.pal.info))
scale.func <- switch(
EXPR = scale.by,
'size' = scale_size,
'radius' = scale_radius,
stop("'scale.by' must be either 'size' or 'radius'")
)
feature.groups <- NULL
if (is.list(features) | any(!is.na(names(features)))) {
feature.groups <- unlist(x = sapply(
X = 1:length(features),
FUN = function(x) {
return(rep(x = names(x = features)[x], each = length(features[[x]])))
}
))
if (any(is.na(x = feature.groups))) {
warning(
"Some feature groups are unnamed.",
call. = FALSE,
immediate. = TRUE
)
}
features <- unlist(x = features)
names(x = feature.groups) <- features
}
cells <- unlist(x = CellsByIdentities(object = object, cells = colnames(object[[assay]]), idents = idents))
data.features <- FetchData(object = object, vars = features, cells = cells)
data.features$id <- if (is.null(x = group.by)) {
Idents(object = object)[cells, drop = TRUE]
} else {
object[[group.by, drop = TRUE]][cells, drop = TRUE]
}
if (!is.factor(x = data.features$id)) {
data.features$id <- factor(x = data.features$id)
}
id.levels <- levels(x = data.features$id)
data.features$id <- as.vector(x = data.features$id)
if (!is.null(x = split.by)) {
splits <- FetchData(object = object, vars = split.by)[cells, split.by]
if (split.colors) {
if (length(x = unique(x = splits)) > length(x = cols)) {
stop(paste0("Need to specify at least ", length(x = unique(x = splits)), " colors using the cols parameter"))
}
cols <- cols[1:length(x = unique(x = splits))]
names(x = cols) <- unique(x = splits)
}
data.features$id <- paste(data.features$id, splits, sep = '_')
unique.splits <- unique(x = splits)
id.levels <- paste0(rep(x = id.levels, each = length(x = unique.splits)), "_", rep(x = unique(x = splits), times = length(x = id.levels)))
}
data.plot <- lapply(
X = unique(x = data.features$id),
FUN = function(ident) {
data.use <- data.features[data.features$id == ident, 1:(ncol(x = data.features) - 1), drop = FALSE]
avg.exp <- apply(
X = data.use,
MARGIN = 2,
FUN = function(x) {
return(mean(x = expm1(x = x)))
}
)
pct.exp <- apply(X = data.use, MARGIN = 2, FUN = PercentAbove, threshold = 0)
return(list(avg.exp = avg.exp, pct.exp = pct.exp))
}
)
names(x = data.plot) <- unique(x = data.features$id)
if (cluster.idents) {
mat <- do.call(
what = rbind,
args = lapply(X = data.plot, FUN = unlist)
)
mat <- scale(x = mat)
id.levels <- id.levels[hclust(d = dist(x = mat))$order]
}
data.plot <- lapply(
X = names(x = data.plot),
FUN = function(x) {
data.use <- as.data.frame(x = data.plot[[x]])
data.use$features.plot <- rownames(x = data.use)
data.use$id <- x
return(data.use)
}
)
data.plot <- do.call(what = 'rbind', args = data.plot)
if (!is.null(x = id.levels)) {
data.plot$id <- factor(x = data.plot$id, levels = id.levels)
}
ngroup <- length(x = levels(x = data.plot$id))
if (ngroup == 1) {
scale <- FALSE
warning(
"Only one identity present, the expression values will be not scaled",
call. = FALSE,
immediate. = TRUE
)
} else if (ngroup < 5 & scale) {
warning(
"Scaling data with a low number of groups may produce misleading results",
call. = FALSE,
immediate. = TRUE
)
}
avg.exp.scaled <- sapply(
X = unique(x = data.plot$features.plot),
FUN = function(x) {
data.use <- data.plot[data.plot$features.plot == x, 'avg.exp']
if (scale) {
data.use <- scale(x = log1p(data.use))
data.use <- MinMax(data = data.use, min = col.min, max = col.max)
} else {
data.use <- log1p(x = data.use)
}
return(data.use)
}
)
avg.exp.scaled <- as.vector(x = t(x = avg.exp.scaled))
if (split.colors) {
avg.exp.scaled <- as.numeric(x = cut(x = avg.exp.scaled, breaks = 20))
}
data.plot$avg.exp.scaled <- avg.exp.scaled
data.plot$features.plot <- factor(
x = data.plot$features.plot,
levels = features
)
data.plot$pct.exp[data.plot$pct.exp < dot.min] <- NA
data.plot$pct.exp <- data.plot$pct.exp * 100
if (split.colors) {
splits.use <- unlist(x = lapply(
X = data.plot$id,
FUN = function(x)
sub(
paste0(".*_(",
paste(sort(unique(x = splits), decreasing = TRUE),
collapse = '|'
),")$"),
"\\1",
x
)
)
)
data.plot$colors <- mapply(
FUN = function(color, value) {
return(colorRampPalette(colors = c('grey', color))(20)[value])
},
color = cols[splits.use],
value = avg.exp.scaled
)
}
color.by <- ifelse(test = split.colors, yes = 'colors', no = 'avg.exp.scaled')
if (!is.na(x = scale.min)) {
data.plot[data.plot$pct.exp < scale.min, 'pct.exp'] <- scale.min
}
if (!is.na(x = scale.max)) {
data.plot[data.plot$pct.exp > scale.max, 'pct.exp'] <- scale.max
}
if (!is.null(x = feature.groups)) {
data.plot$feature.groups <- factor(
x = feature.groups[data.plot$features.plot],
levels = unique(x = feature.groups)
)
}
plot <- ggplot(data = data.plot, mapping = aes(x = .data[["features.plot"]], y = .data[["id"]])) +
geom_point(mapping = aes(size = .data[["pct.exp"]], color = .data[[color.by]])) +
scale.func(range = c(0, dot.scale), limits = c(scale.min, scale.max)) +
theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
guides(size = guide_legend(title = 'Percent Expressed')) +
labs(
x = 'Features',
y = ifelse(test = is.null(x = split.by), yes = 'Identity', no = 'Split Identity')
) +
theme_cowplot()
if (!is.null(x = feature.groups)) {
plot <- plot + facet_grid(
rows = ~feature.groups,
scales = "free_x",
space = "free_x",
switch = "y"
) + theme(
panel.spacing = unit(x = 1, units = "lines"),
strip.background = element_blank()
)
}
if (split.colors) {
plot <- plot + scale_color_identity()
} else if (length(x = cols) == 1) {
plot <- plot + scale_color_distiller(palette = cols)
} else {
plot <- plot + scale_color_gradient(low = cols[1], high = cols[2])
}
if (!split.colors) {
plot <- plot + guides(color = guide_colorbar(title = 'Average Expression'))
}
return(plot)
}
#' Quickly Pick Relevant Dimensions
#'
#' Plots per-component standard deviations (or approximate singular values if running PCAFast),
#' percent variance explained per principal component, or cumulative percent variance explained,
#' to help pick an elbow in the graph. This elbow often corresponds well with significant
#' dimensions and is much faster to run than Jackstraw.
#'
#' @param object Seurat object
#' @param ndims Number of dimensions to plot (positive integer; capped by stored components)
#' @param reduction Reduction technique to plot (default is 'pca')
#' @param plot_type One of \code{"stdev"} (default), \code{"variance"} (per-PC \% variance), or
#' \code{"cumulative_variance"} (running sum of those percentages; equals 100\% at the last
#' stored PC when \code{ndims} spans all of them)
#'
#' @return A ggplot object
#'
#' @importFrom cowplot theme_cowplot
#' @importFrom ggplot2 ggplot geom_point labs aes
#' @export
#' @concept visualization
#'
#' @examples
#' data("pbmc_small")
#' ElbowPlot(object = pbmc_small)
#' ElbowPlot(object = pbmc_small, plot_type = "variance")
#' ElbowPlot(object = pbmc_small, plot_type = "cumulative_variance")
#'
ElbowPlot <- function(object, ndims = 20, reduction = 'pca', plot_type = c("stdev", "variance", "cumulative_variance")) {
plot_type <- match.arg(plot_type)
if (!is.numeric(ndims) || length(ndims) != 1L || !is.finite(ndims) || ndims < 1 || ndims != as.integer(ndims)) {
stop("'ndims' must be a single positive integer", call. = FALSE)
}
ndims <- as.integer(ndims)
data.use <- Stdev(object = object, reduction = reduction)
if (length(x = data.use) == 0) {
stop(paste("No standard deviation info stored for", reduction))
}
if (anyNA(data.use)) {
stop("Standard deviations contain NA for reduction ", reduction, call. = FALSE)
}
if (ndims > length(x = data.use)) {
warning("The object only has information for ", length(x = data.use), " dimensions")
ndims <- length(x = data.use)
}
if (plot_type == "stdev") {
y_label <- "Standard Deviation"
y_data <- data.use[1:ndims]
} else {
den <- sum(data.use^2)
if (!is.finite(den) || den == 0) {
stop("Cannot compute variance explained: sum of squared standard deviations is not positive for reduction ", reduction, call. = FALSE)
}
pct <- data.use^2 / den * 100
if (plot_type == "variance") {
y_label <- "Percentage of Variance Explained"
y_data <- pct[1:ndims]
} else {
y_label <- "Cumulative % Variance Explained"
y_data <- cumsum(pct)[1:ndims]
}
}
plot <- ggplot(data = data.frame(dims = 1:ndims, y_data = y_data)) +
geom_point(mapping = aes(x = .data[["dims"]], y = .data[["y_data"]])) +
labs(x = gsub(pattern = '_$', replacement = '', x = Key(object = object[[reduction]])), y = y_label) +
theme_cowplot()
return(plot)
}
#' Boxplot of correlation of a variable (e.g. number of UMIs) with expression
#' data
#'
#' @param object Seurat object
#' @param assay Assay where the feature grouping info and correlations are
#' stored
#' @param feature.group Name of the column in meta.features where the feature
#' grouping info is stored
#' @param cor Name of the column in meta.features where correlation info is
#' stored
#'
#' @return Returns a ggplot boxplot of correlations split by group
#'
#' @importFrom ggplot2 geom_boxplot scale_fill_manual geom_hline
#' @importFrom cowplot theme_cowplot
#' @importFrom scales brewer_pal
#' @importFrom stats complete.cases
#'
#' @export
#' @concept visualization
#'
GroupCorrelationPlot <- function(
object,
assay = NULL,
feature.group = "feature.grp",
cor = "nCount_RNA_cor"
) {
assay <- assay %||% DefaultAssay(object = object)
data <- object[[assay]][[c(feature.group, cor)]]
data <- data[complete.cases(data), ]
colnames(x = data) <- c('grp', 'cor')
data$grp <- as.character(data$grp)
plot <- ggplot(data = data, aes(x = .data[["grp"]], y = .data[["cor"]], fill = .data[["grp"]])) +
geom_boxplot() +
theme_cowplot() +
scale_fill_manual(values = rev(x = brewer_pal(palette = 'YlOrRd')(n = 7))) +
ylab(paste(
"Correlation with",
gsub(x = cor, pattern = "_cor", replacement = "")
)) +
geom_hline(yintercept = 0) +
NoLegend() +
theme(
axis.line.x = element_blank(),
axis.title.x = element_blank(),
axis.ticks.x = element_blank(),
axis.text.x = element_blank()
)
return(plot)
}
#' JackStraw Plot
#'
#' Plots the results of the JackStraw analysis for PCA significance. For each
#' PC, plots a QQ-plot comparing the distribution of p-values for all genes
#' across each PC, compared with a uniform distribution. Also determines a
#' p-value for the overall significance of each PC (see Details).
#'
#' Significant PCs should show a p-value distribution (black curve) that is
#' strongly skewed to the left compared to the null distribution (dashed line)
#' The p-value for each PC is based on a proportion test comparing the number
#' of genes with a p-value below a particular threshold (score.thresh), compared with the
#' proportion of genes expected under a uniform distribution of p-values.
#'
#' @param object Seurat object
#' @param dims Dims to plot
#' @param cols Vector of colors, each color corresponds to an individual PC. This may also be a single character
#' or numeric value corresponding to a palette as specified by \code{\link[RColorBrewer]{brewer.pal.info}}.
#' By default, ggplot2 assigns colors. We also include a number of palettes from the pals package.
#' See \code{\link{DiscretePalette}} for details.
#' @param reduction reduction to pull jackstraw info from
#' @param xmax X-axis maximum on each QQ plot.
#' @param ymax Y-axis maximum on each QQ plot.
#'
#' @return A ggplot object
#'
#' @author Omri Wurtzel
#' @seealso \code{\link{ScoreJackStraw}}
#'
#' @importFrom stats qunif
#' @importFrom scales hue_pal
#' @importFrom ggplot2 ggplot stat_qq labs xlim ylim
#' coord_flip geom_abline guides guide_legend
#' @importFrom cowplot theme_cowplot
#'
#' @export
#' @concept visualization
#'
#' @examples
#' data("pbmc_small")
#' JackStrawPlot(object = pbmc_small)
#'
JackStrawPlot <- function(
object,
dims = 1:5,
cols = NULL,
reduction = 'pca',
xmax = 0.1,
ymax = 0.3
) {
pAll <- JS(object = object[[reduction]], slot = 'empirical')
if (max(dims) > ncol(x = pAll)) {
stop("Max dimension is ", ncol(x = pAll))
}
pAll <- pAll[, dims, drop = FALSE]
pAll <- as.data.frame(x = pAll)
data.plot <- Melt(x = pAll)
colnames(x = data.plot) <- c("Contig", "PC", "Value")
score.df <- JS(object = object[[reduction]], slot = 'overall')
if (nrow(x = score.df) < max(dims)) {
stop("Jackstraw procedure not scored for all the provided dims. Please run ScoreJackStraw.")
}
score.df <- score.df[dims, , drop = FALSE]
if (nrow(x = score.df) == 0) {
stop(paste0("JackStraw hasn't been scored. Please run ScoreJackStraw before plotting."))
}
data.plot$PC.Score <- rep(
x = paste0("PC ", score.df[ ,"PC"], ": ", sprintf("%1.3g", score.df[ ,"Score"])),
each = length(x = unique(x = data.plot$Contig))
)
data.plot$PC.Score <- factor(
x = data.plot$PC.Score,
levels = paste0("PC ", score.df[, "PC"], ": ", sprintf("%1.3g", score.df[, "Score"]))
)
if (is.null(x = cols)) {
cols <- hue_pal()(length(x = dims))
}
if (length(x = cols) < length(x = dims)) {
stop("Not enough colors for the number of dims selected")
}
gp <- ggplot(data = data.plot, mapping = aes(sample = .data[['Value']], color = .data[['PC.Score']])) +
stat_qq(distribution = qunif) +
labs(x = "Theoretical [runif(1000)]", y = "Empirical") +
scale_color_manual(values = cols) +
xlim(0, ymax) +
ylim(0, xmax) +
coord_flip() +
geom_abline(intercept = 0, slope = 1, linetype = "dashed", na.rm = TRUE) +
guides(color = guide_legend(title = "PC: p-value")) +
theme_cowplot()
return(gp)
}
#' Plot clusters as a tree
#'
#' Plots previously computed tree (from BuildClusterTree)
#'
#' @param object Seurat object
#' @param direction A character string specifying the direction of the tree (default is downwards)
#' Possible options: "rightwards", "leftwards", "upwards", and "downwards".
#' @param \dots Additional arguments to
#' \code{\link[ape:plot.phylo]{ape::plot.phylo}}
#'
#' @return Plots dendogram (must be precomputed using BuildClusterTree), returns no value
#'
#' @export
#' @concept visualization
#'
#' @examples
#' \dontrun{
#' if (requireNamespace("ape", quietly = TRUE)) {
#' data("pbmc_small")
#' pbmc_small <- BuildClusterTree(object = pbmc_small)
#' PlotClusterTree(object = pbmc_small)
#' }
#' }
PlotClusterTree <- function(object, direction = "downwards", ...) {
if (isFALSE(x = requireNamespace('ape', quietly = TRUE))) {
stop(cluster.ape, call. = FALSE)
}
if (is.null(x = Tool(object = object, slot = "BuildClusterTree"))) {
stop("Phylogenetic tree does not exist, build using BuildClusterTree")
}
data.tree <- Tool(object = object, slot = "BuildClusterTree")
ape::plot.phylo(x = data.tree, direction = direction, ...)
ape::nodelabels()
}
#' Visualize Dimensional Reduction genes
#'
#' Visualize top genes associated with reduction components
#'
#' @param object Seurat object
#' @param reduction Reduction technique to visualize results for
#' @param dims Number of dimensions to display
#' @param nfeatures Number of genes to display
#' @param col Color of points to use
#' @param projected Use reduction values for full dataset (i.e. projected
#' dimensional reduction values)
#' @param balanced Return an equal number of genes with + and - scores. If
#' FALSE (default), returns the top genes ranked by the scores absolute values
#' @param ncol Number of columns to display
#' @param combine Combine plots into a single \code{patchwork}
#' ggplot object. If \code{FALSE}, return a list of ggplot objects
#'
#' @return A \code{patchwork} ggplot object if
#' \code{combine = TRUE}; otherwise, a list of ggplot objects
#'
#' @importFrom patchwork wrap_plots
#' @importFrom cowplot theme_cowplot
#' @importFrom ggplot2 ggplot geom_point labs
#' @export
#' @concept visualization
#'
#' @examples
#' data("pbmc_small")
#' VizDimLoadings(object = pbmc_small)
#'
VizDimLoadings <- function(
object,
dims = 1:5,
nfeatures = 30,
col = 'blue',
reduction = 'pca',
projected = FALSE,
balanced = FALSE,
ncol = NULL,
combine = TRUE
) {
if (is.null(x = ncol)) {
ncol <- 2
if (length(x = dims) == 1) {
ncol <- 1
}
if (length(x = dims) > 6) {
ncol <- 3
}
if (length(x = dims) > 9) {
ncol <- 4
}
}
loadings <- Loadings(object = object[[reduction]], projected = projected)
features <- lapply(
X = dims,
FUN = TopFeatures,
object = object[[reduction]],
nfeatures = nfeatures,
projected = projected,
balanced = balanced
)
features <- lapply(
X = features,
FUN = unlist,
use.names = FALSE
)
loadings <- loadings[unlist(x = features), dims, drop = FALSE]
names(x = features) <- colnames(x = loadings) <- as.character(x = dims)
plots <- lapply(
X = as.character(x = dims),
FUN = function(i) {
data.plot <- as.data.frame(x = loadings[features[[i]], i, drop = FALSE])
colnames(x = data.plot) <- paste0(Key(object = object[[reduction]]), i)
data.plot$feature <- factor(x = rownames(x = data.plot), levels = rownames(x = data.plot))
plot <- ggplot(
data = data.plot,
mapping = aes(x = .data[[paste0(Key(object = object[[reduction]]), i)]], y = .data[['feature']])
) +
geom_point(col = col) +
labs(y = NULL) + theme_cowplot()
return(plot)
}
)
if (combine) {
plots <- wrap_plots(plots, ncol = ncol)
}
return(plots)
}
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Exported utility functions
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#' Augments ggplot2-based plot with a PNG image.
#'
#' Creates "vector-friendly" plots. Does this by saving a copy of the plot as a PNG file,
#' then adding the PNG image with \code{\link[ggplot2]{annotation_raster}} to a blank plot
#' of the same dimensions as \code{plot}. Please note: original legends and axes will be lost
#' during augmentation.
#'
#' @param plot A ggplot object
#' @param width,height Width and height of PNG version of plot
#' @param dpi Plot resolution
#'
#' @return A ggplot object
#'
#' @importFrom png readPNG
#' @importFrom ggplot2 ggplot_build ggsave ggplot geom_blank annotation_raster ggtitle
#'
#' @export
#' @concept visualization
#'
#' @examples
#' \dontrun{
#' data("pbmc_small")
#' plot <- DimPlot(object = pbmc_small)
#' AugmentPlot(plot = plot)
#' }
#'
AugmentPlot <- function(plot, width = 10, height = 10, dpi = 100) {
pbuild.params <- ggplot_build(plot = plot)$layout$panel_params[[1]]
range.values <- c(
pbuild.params$x.range,
pbuild.params$y.range
)
xyparams <- GetXYAesthetics(
plot = plot,
geom = class(x = plot$layers[[1]]$geom)[1]
)
title <- plot$labels$title
tmpfile <- tempfile(fileext = '.png')
ggsave(
filename = tmpfile,
plot = plot + NoLegend() + NoAxes() + theme(plot.title = element_blank()),
width = width,
height = height,
dpi = dpi
)
img <- readPNG(source = tmpfile)
file.remove(tmpfile)
blank <- ggplot(
data = plot$data,
mapping = aes(x = .data[[xyparams$x]], y = .data[[xyparams$y]])
) + geom_blank()
blank <- blank + plot$theme + ggtitle(label = title)
blank <- blank + annotation_raster(
raster = img,
xmin = range.values[1],
xmax = range.values[2],
ymin = range.values[3],
ymax = range.values[4]
)
return(blank)
}
#' Automagically calculate a point size for ggplot2-based scatter plots
#'
#' It happens to look good
#'
#' @param data A data frame being passed to ggplot2
#' @param raster If TRUE, point size is set to 1
#'
#' @return The "optimal" point size for visualizing these data
#'
#' @export
#' @concept visualization
#'
#' @examples
#' df <- data.frame(x = rnorm(n = 10000), y = runif(n = 10000))
#' AutoPointSize(data = df)
#'
AutoPointSize <- function(data, raster = NULL) {
return(ifelse(
test = isTRUE(x = raster),
yes = 1,
no = min(1583 / nrow(x = data), 1)
))
}
#' Determine text color based on background color
#'
#' @param background A vector of background colors; supports R color names and
#' hexadecimal codes
#' @param threshold Intensity threshold for light/dark cutoff; intensities
#' greater than \code{theshold} yield \code{dark}, others yield \code{light}
#' @param w3c Use \href{https://www.w3.org/TR/WCAG20/}{W3C} formula for calculating
#' background text color; ignores \code{threshold}
#' @param dark Color for dark text
#' @param light Color for light text
#'
#' @return A named vector of either \code{dark} or \code{light}, depending on
#' \code{background}; names of vector are \code{background}
#'
#' @export
#' @concept visualization
#'
#' @source \url{https://stackoverflow.com/questions/3942878/how-to-decide-font-color-in-white-or-black-depending-on-background-color}
#'
#' @examples
#' BGTextColor(background = c('black', 'white', '#E76BF3'))
#'
BGTextColor <- function(
background,
threshold = 186,
w3c = FALSE,
dark = 'black',
light = 'white'
) {
if (w3c) {
luminance <- Luminance(color = background)
threshold <- 179
return(ifelse(
test = luminance > sqrt(x = 1.05 * 0.05) - 0.05,
yes = dark,
no = light
))
}
return(ifelse(
test = Intensity(color = background) > threshold,
yes = dark,
no = light
))
}
#' @export
#' @concept visualization
#'
#' @rdname CustomPalette
#' @aliases BlackAndWhite
#'
#' @examples
#' df <- data.frame(x = rnorm(n = 100, mean = 20, sd = 2), y = rbinom(n = 100, size = 100, prob = 0.2))
#' plot(df, col = BlackAndWhite())
#'
BlackAndWhite <- function(mid = NULL, k = 50) {
return(CustomPalette(low = "white", high = "black", mid = mid, k = k))
}
#' @export
#' @concept visualization
#'
#' @rdname CustomPalette
#' @aliases BlueAndRed
#'
#' @examples
#' df <- data.frame(x = rnorm(n = 100, mean = 20, sd = 2), y = rbinom(n = 100, size = 100, prob = 0.2))
#' plot(df, col = BlueAndRed())
#'
BlueAndRed <- function(k = 50) {
return(CustomPalette(low = "#313695" , high = "#A50026", mid = "#FFFFBF", k = k))
}
#' Cell Selector
#'
#' Select points on a scatterplot and get information about them
#'
#' @param plot A ggplot2 plot
#' @param object An optional Seurat object; if passes, will return an object
#' with the identities of selected cells set to \code{ident}
#' @param ident An optional new identity class to assign the selected cells
#' @param ... Ignored
#'
#' @return If \code{object} is \code{NULL}, the names of the points selected;
#' otherwise, a Seurat object with the selected cells identity classes set to
#' \code{ident}
#'
#' @importFrom miniUI miniPage gadgetTitleBar miniTitleBarButton
#' miniContentPanel
#' @importFrom shiny fillRow plotOutput brushOpts reactiveValues observeEvent
#' stopApp brushedPoints renderPlot runGadget
#'
#' @export
#' @concept visualization
#'
#' @seealso \code{\link{DimPlot}} \code{\link{FeaturePlot}}
#'
#' @examples
#' \dontrun{
#' data("pbmc_small")
#' plot <- DimPlot(object = pbmc_small)
#' # Follow instructions in the terminal to select points
#' cells.located <- CellSelector(plot = plot)
#' cells.located
#' # Automatically set the identity class of selected cells and return a new Seurat object
#' pbmc_small <- CellSelector(plot = plot, object = pbmc_small, ident = 'SelectedCells')
#' }
#'
CellSelector <- function(plot, object = NULL, ident = 'SelectedCells', ...) {
# Set up the gadget UI
ui <- miniPage(
gadgetTitleBar(
title = "Cell Selector",
left = miniTitleBarButton(inputId = "reset", label = "Reset")
),
miniContentPanel(
fillRow(
plotOutput(
outputId = "plot",
height = '100%',
brush = brushOpts(
id = 'brush',
delay = 100,
delayType = 'debounce',
clip = TRUE,
resetOnNew = FALSE
)
)
),
)
)
# Get some plot information
if (inherits(x = plot, what = 'patchwork')) {
if (length(x = plot$patches$plots)) {
warning(
"Multiple plots passed, using last plot",
call. = FALSE,
immediate. = TRUE
)
}
class(x = plot) <- grep(
pattern = 'patchwork',
x = class(x = plot),
value = TRUE,
invert = TRUE
)
}
xy.aes <- GetXYAesthetics(plot = plot)
dark.theme <- !is.null(x = plot$theme$plot.background$fill) &&
plot$theme$plot.background$fill == 'black'
plot.data <- GGpointToBase(plot = plot, do.plot = FALSE)
plot.data$selected_ <- FALSE
rownames(x = plot.data) <- rownames(x = plot$data)
colnames(x = plot.data) <- gsub(
pattern = '-',
replacement = '.',
x = colnames(x = plot.data)
)
# Server function
server <- function(input, output, session) {
plot.env <- reactiveValues(data = plot.data)
# Event handlers
observeEvent(
eventExpr = input$done,
handlerExpr = {
PlotBuild(data = plot.env$data, dark.theme = dark.theme)
selected <- rownames(x = plot.data)[plot.env$data$selected_]
if (inherits(x = object, what = 'Seurat')) {
if (!all(selected %in% Cells(x = object))) {
stop("Cannot find the selected cells in the Seurat object, please be sure you pass the same object used to generate the plot")
}
Idents(object = object, cells = selected) <- ident
selected <- object
}
stopApp(returnValue = selected)
}
)
observeEvent(
eventExpr = input$reset,
handlerExpr = {
plot.env$data <- plot.data
session$resetBrush(brushId = 'brush')
}
)
observeEvent(
eventExpr = input$brush,
handlerExpr = {
plot.env$data <- brushedPoints(
df = plot.data,
brush = input$brush,
xvar = xy.aes$x,
yvar = xy.aes$y,
allRows = TRUE
)
plot.env$data$color <- ifelse(
test = plot.env$data$selected_,
yes = '#DE2D26',
no = '#C3C3C3'
)
}
)
# Render the plot
output$plot <- renderPlot(expr = PlotBuild(
data = plot.env$data,
dark.theme = dark.theme
))
}
return(runGadget(app = ui, server = server))
}
#' Move outliers towards center on dimension reduction plot
#'
#' @param object Seurat object
#' @param reduction Name of DimReduc to adjust
#' @param dims Dimensions to visualize
#' @param group.by Group (color) cells in different ways (for example, orig.ident)
#' @param outlier.sd Controls the outlier distance
#' @param reduction.key Key for DimReduc that is returned
#'
#' @return Returns a DimReduc object with the modified embeddings
#'
#' @export
#' @concept visualization
#'
#' @examples
#' \dontrun{
#' data("pbmc_small")
#' pbmc_small <- FindClusters(pbmc_small, resolution = 1.1)
#' pbmc_small <- RunUMAP(pbmc_small, dims = 1:5)
#' DimPlot(pbmc_small, reduction = "umap")
#' pbmc_small[["umap_new"]] <- CollapseEmbeddingOutliers(pbmc_small,
#' reduction = "umap", reduction.key = 'umap_', outlier.sd = 0.5)
#' DimPlot(pbmc_small, reduction = "umap_new")
#' }
#'
CollapseEmbeddingOutliers <- function(
object,
reduction = 'umap',
dims = 1:2,
group.by = 'ident',
outlier.sd = 2,
reduction.key = 'UMAP_'
) {
embeddings <- Embeddings(object = object[[reduction]])[, dims]
idents <- FetchData(object = object, vars = group.by)
data.medians <- sapply(X = dims, FUN = function(x) {
tapply(X = embeddings[, x], INDEX = idents, FUN = median)
})
data.sd <- apply(X = data.medians, MARGIN = 2, FUN = sd)
data.medians.scale <- as.matrix(x = scale(x = data.medians, center = TRUE, scale = TRUE))
data.medians.scale[abs(x = data.medians.scale) < outlier.sd] <- 0
data.medians.scale <- sign(x = data.medians.scale) * (abs(x = data.medians.scale) - outlier.sd)
data.correct <- Sweep(
x = data.medians.scale,
MARGIN = 2,
STATS = data.sd,
FUN = "*"
)
data.correct <- data.correct[abs(x = apply(X = data.correct, MARGIN = 1, FUN = min)) > 0, ]
new.embeddings <- embeddings
for (i in rownames(x = data.correct)) {
cells.correct <- rownames(x = idents)[idents[, "ident"] == i]
new.embeddings[cells.correct, ] <- Sweep(
x = new.embeddings[cells.correct,],
MARGIN = 2,
STATS = data.correct[i, ],
FUN = "-"
)
}
reduc <- CreateDimReducObject(
embeddings = new.embeddings,
loadings = Loadings(object = object[[reduction]]),
assay = slot(object = object[[reduction]], name = "assay.used"),
key = reduction.key
)
return(reduc)
}
#' Combine ggplot2-based plots into a single plot
#'
#' @param plots A list of gg objects
#' @param ncol Number of columns
#' @param legend Combine legends into a single legend
#' choose from 'right' or 'bottom'; pass 'none' to remove legends, or \code{NULL}
#' to leave legends as they are
#' @param ... Extra parameters passed to plot_grid
#'
#' @return A combined plot
#'
#' @importFrom cowplot plot_grid get_legend
#' @export
#' @concept visualization
#'
#' @examples
#' data("pbmc_small")
#' pbmc_small[['group']] <- sample(
#' x = c('g1', 'g2'),
#' size = ncol(x = pbmc_small),
#' replace = TRUE
#' )
#' plot1 <- FeaturePlot(
#' object = pbmc_small,
#' features = 'MS4A1',
#' split.by = 'group'
#' )
#' plot2 <- FeaturePlot(
#' object = pbmc_small,
#' features = 'FCN1',
#' split.by = 'group'
#' )
#' CombinePlots(
#' plots = list(plot1, plot2),
#' legend = 'none',
#' nrow = length(x = unique(x = pbmc_small[['group', drop = TRUE]]))
#' )
#'
CombinePlots <- function(plots, ncol = NULL, legend = NULL, ...) {
.Deprecated(msg = "CombinePlots is being deprecated. Plots should now be combined using the patchwork system.")
plots.combined <- if (length(x = plots) > 1) {
if (!is.null(x = legend)) {
if (legend != 'none') {
plot.legend <- get_legend(plot = plots[[1]] + theme(legend.position = legend))
}
plots <- lapply(
X = plots,
FUN = function(x) {
return(x + NoLegend())
}
)
}
plots.combined <- plot_grid(
plotlist = plots,
ncol = ncol,
align = 'hv',
...
)
if (!is.null(x = legend)) {
plots.combined <- switch(
EXPR = legend,
'bottom' = plot_grid(
plots.combined,
plot.legend,
ncol = 1,
rel_heights = c(1, 0.2)
),
'right' = plot_grid(
plots.combined,
plot.legend,
rel_widths = c(3, 0.3)
),
plots.combined
)
}
plots.combined
} else {
plots[[1]]
}
return(plots.combined)
}
#' Create a custom color palette
#'
#' Creates a custom color palette based on low, middle, and high color values
#'
#' @param low low color
#' @param high high color
#' @param mid middle color. Optional.
#' @param k number of steps (colors levels) to include between low and high values
#'
#' @return A color palette for plotting
#'
#' @importFrom grDevices col2rgb rgb
#' @export
#' @concept visualization
#'
#' @rdname CustomPalette
#' @examples
#' myPalette <- CustomPalette()
#' myPalette
#'
CustomPalette <- function(
low = "white",
high = "red",
mid = NULL,
k = 50
) {
low <- col2rgb(col = low) / 255
high <- col2rgb(col = high) / 255
if (is.null(x = mid)) {
r <- seq(from = low[1], to = high[1], len = k)
g <- seq(from = low[2], to = high[2], len = k)
b <- seq(from = low[3], to = high[3], len = k)
} else {
k2 <- round(x = k / 2)
mid <- col2rgb(col = mid) / 255
r <- c(
seq(from = low[1], to = mid[1], len = k2),
seq(from = mid[1], to = high[1], len = k2)
)
g <- c(
seq(from = low[2], to = mid[2], len = k2),
seq(from = mid[2], to = high[2],len = k2)
)
b <- c(
seq(from = low[3], to = mid[3], len = k2),
seq(from = mid[3], to = high[3], len = k2)
)
}
return(rgb(red = r, green = g, blue = b))
}
#' Discrete colour palettes from pals
#'
#' These are included here because pals depends on a number of compiled
#' packages, and this can lead to increases in run time for Travis,
#' and generally should be avoided when possible.
#'
#' These palettes are a much better default for data with many classes
#' than the default ggplot2 palette.
#'
#' Many thanks to Kevin Wright for writing the pals package.
#'
#' @param n Number of colours to be generated.
#' @param palette Options are
#' "alphabet", "alphabet2", "glasbey", "polychrome", "stepped", and "parade".
#' Can be omitted and the function will use the one based on the requested n.
#' @param shuffle Shuffle the colors in the selected palette.
#'
#' @return A vector of colors
#'
#' @details
#' Taken from the pals package (Licence: GPL-3).
#' \url{https://cran.r-project.org/package=pals}
#' Credit: Kevin Wright
#'
#' @export
#' @concept visualization
#'
DiscretePalette <- function(n, palette = NULL, shuffle = FALSE) {
palettes <- list(
alphabet = c(
"#F0A0FF", "#0075DC", "#993F00", "#4C005C", "#191919", "#005C31",
"#2BCE48", "#FFCC99", "#808080", "#94FFB5", "#8F7C00", "#9DCC00",
"#C20088", "#003380", "#FFA405", "#FFA8BB", "#426600", "#FF0010",
"#5EF1F2", "#00998F", "#E0FF66", "#740AFF", "#990000", "#FFFF80",
"#FFE100", "#FF5005"
),
alphabet2 = c(
"#AA0DFE", "#3283FE", "#85660D", "#782AB6", "#565656", "#1C8356",
"#16FF32", "#F7E1A0", "#E2E2E2", "#1CBE4F", "#C4451C", "#DEA0FD",
"#FE00FA", "#325A9B", "#FEAF16", "#F8A19F", "#90AD1C", "#F6222E",
"#1CFFCE", "#2ED9FF", "#B10DA1", "#C075A6", "#FC1CBF", "#B00068",
"#FBE426", "#FA0087"
),
glasbey = c(
"#0000FF", "#FF0000", "#00FF00", "#000033", "#FF00B6", "#005300",
"#FFD300", "#009FFF", "#9A4D42", "#00FFBE", "#783FC1", "#1F9698",
"#FFACFD", "#B1CC71", "#F1085C", "#FE8F42", "#DD00FF", "#201A01",
"#720055", "#766C95", "#02AD24", "#C8FF00", "#886C00", "#FFB79F",
"#858567", "#A10300", "#14F9FF", "#00479E", "#DC5E93", "#93D4FF",
"#004CFF", "#F2F318"
),
polychrome = c(
"#5A5156", "#E4E1E3", "#F6222E", "#FE00FA", "#16FF32", "#3283FE",
"#FEAF16", "#B00068", "#1CFFCE", "#90AD1C", "#2ED9FF", "#DEA0FD",
"#AA0DFE", "#F8A19F", "#325A9B", "#C4451C", "#1C8356", "#85660D",
"#B10DA1", "#FBE426", "#1CBE4F", "#FA0087", "#FC1CBF", "#F7E1A0",
"#C075A6", "#782AB6", "#AAF400", "#BDCDFF", "#822E1C", "#B5EFB5",
"#7ED7D1", "#1C7F93", "#D85FF7", "#683B79", "#66B0FF", "#3B00FB"
),
stepped = c(
"#990F26", "#B33E52", "#CC7A88", "#E6B8BF", "#99600F", "#B3823E",
"#CCAA7A", "#E6D2B8", "#54990F", "#78B33E", "#A3CC7A", "#CFE6B8",
"#0F8299", "#3E9FB3", "#7ABECC", "#B8DEE6", "#3D0F99", "#653EB3",
"#967ACC", "#C7B8E6", "#333333", "#666666", "#999999", "#CCCCCC"
),
parade = c(
'#ff6969', '#9b37ff', '#cd3737', '#69cdff', '#ffff69', '#69cdcd',
'#9b379b', '#3737cd', '#ffff9b', '#cdff69', '#ff9b37', '#37ffff',
'#9b69ff', '#37cd69', '#ff3769', '#ff3737', '#37ff9b', '#cdcd37',
'#3769cd', '#37cdff', '#9b3737', '#ff699b', '#9b9bff', '#cd9b37',
'#69ff37', '#cd3769', '#cd69cd', '#cd6937', '#3737ff', '#cdcd69',
'#ff9b69', '#cd37cd', '#9bff37', '#cd379b', '#cd6969', '#69ff9b',
'#ff379b', '#9bff9b', '#6937ff', '#69cd37', '#cdff37', '#9bff69',
'#9b37cd', '#ff37ff', '#ff37cd', '#ffff37', '#37cd9b', '#379bff',
'#ffcd37', '#379b37', '#ff9bff', '#379b9b', '#69ffcd', '#379bcd',
'#ff69ff', '#ff9b9b', '#37ff69', '#ff6937', '#6969ff', '#699bff',
'#ffcd69', '#69ffff', '#37ff37', '#6937cd', '#37cd37', '#3769ff',
'#cd69ff', '#6969cd', '#9bcd37', '#69ff69', '#37cdcd', '#cd37ff',
'#37379b', '#37ffcd', '#69cd69', '#ff69cd', '#9bffff', '#9b9b37'
)
)
if (is.null(x = n)) {
return(names(x = palettes))
}
if (is.null(x = palette)) {
if (n <= 26) {
palette <- "alphabet"
} else if (n <= 32) {
palette <- "glasbey"
} else {
palette <- "polychrome"
}
}
palette.vec <- palettes[[palette]]
if (n > length(x = palette.vec)) {
warning("Not enough colours in specified palette")
}
if (isTRUE(shuffle)) {
palette.vec <- sample(palette.vec)
}
palette <- palette.vec[seq_len(length.out = n)]
return(palette)
}
#' @rdname CellSelector
#' @export
#' @concept visualization
#'
FeatureLocator <- function(plot, ...) {
.Defunct(
new = 'CellSelector',
package = 'Seurat',
msg = "'FeatureLocator' has been replaced by 'CellSelector'"
)
}
#' Hover Locator
#'
#' Get quick information from a scatterplot by hovering over points
#'
#' @param plot A ggplot2 plot
#' @param information An optional dataframe or matrix of extra information to be displayed on hover
#' @param dark.theme Plot using a dark theme?
#' @param axes Display or hide x- and y-axes
#' @param ... Extra parameters to be passed to \code{\link[plotly]{layout}}
#'
#' @importFrom ggplot2 ggplot_build
#' @importFrom plotly plot_ly layout add_annotations
#' @export
#' @concept visualization
#'
#' @seealso \code{\link[plotly]{layout}} \code{\link[ggplot2]{ggplot_build}}
#' \code{\link{DimPlot}} \code{\link{FeaturePlot}}
#'
#' @examples
#' \dontrun{
#' data("pbmc_small")
#' plot <- DimPlot(object = pbmc_small)
#' HoverLocator(plot = plot, information = FetchData(object = pbmc_small, vars = 'percent.mito'))
#' }
#'
HoverLocator <- function(
plot,
information = NULL,
axes = TRUE,
dark.theme = FALSE,
...
) {
# Use GGpointToBase because we already have ggplot objects
# with colors (which are annoying in plotly)
plot.build <- suppressWarnings(expr = GGpointToPlotlyBuild(
plot = plot,
information = information,
...
))
data <- ggplot_build(plot = plot)$plot$data
# Set up axis labels here
# Also, a bunch of stuff to get axis lines done properly
if (axes) {
xaxis <- list(
title = names(x = data)[1],
showgrid = FALSE,
zeroline = FALSE,
showline = TRUE
)
yaxis <- list(
title = names(x = data)[2],
showgrid = FALSE,
zeroline = FALSE,
showline = TRUE
)
} else {
xaxis <- yaxis <- list(visible = FALSE)
}
# Check for dark theme
if (dark.theme) {
title <- list(color = 'white')
xaxis <- c(xaxis, color = 'white')
yaxis <- c(yaxis, color = 'white')
plotbg <- 'black'
} else {
title = list(color = 'black')
plotbg = 'white'
}
# The `~' means pull from the data passed (this is why we reset the names)
# Use I() to get plotly to accept the colors from the data as is
# Set hoverinfo to 'text' to override the default hover information
# rather than append to it
p <- layout(
p = plot_ly(
data = plot.build,
x = ~x,
y = ~y,
type = 'scatter',
mode = 'markers',
color = ~I(color),
hoverinfo = 'text',
text = ~feature
),
xaxis = xaxis,
yaxis = yaxis,
title = plot$labels$title,
titlefont = title,
paper_bgcolor = plotbg,
plot_bgcolor = plotbg,
...
)
# Add labels
label.layer <- which(x = sapply(
X = plot$layers,
FUN = function(x) {
return(inherits(x = x$geom, what = c('GeomText', 'GeomTextRepel')))
}
))
if (length(x = label.layer) == 1) {
p <- add_annotations(
p = p,
x = plot$layers[[label.layer]]$data[, 1],
y = plot$layers[[label.layer]]$data[, 2],
xref = "x",
yref = "y",
text = plot$layers[[label.layer]]$data[, 3],
xanchor = 'right',
showarrow = FALSE,
font = list(size = plot$layers[[label.layer]]$aes_params$size * 4)
)
}
return(p)
}
#' Get the intensity and/or luminance of a color
#'
#' @param color A vector of colors
#'
#' @return A vector of intensities/luminances for each color
#'
#' @name contrast-theory
#' @rdname contrast-theory
#'
#' @importFrom grDevices col2rgb
#'
#' @export
#' @concept visualization
#'
#' @source \url{https://stackoverflow.com/questions/3942878/how-to-decide-font-color-in-white-or-black-depending-on-background-color}
#'
#' @examples
#' Intensity(color = c('black', 'white', '#E76BF3'))
#'
Intensity <- function(color) {
intensities <- apply(
X = col2rgb(col = color),
MARGIN = 2,
FUN

visualization.R at commit 586015a, under other · at the source

Overview

Authors: Dou Ye^1,2, Haipeng Zhou², Suqing Qu², Zhaoyan Wang², Fan Zhang², Xiaohua Wang², Yuan Zhao², Jialan Liang², Qian Wang², Zuo Luan^2,3, Yinxiang Yang²

ORCID iDs: Dou Ye, Haipeng Zhou, Suqing Qu, Zhaoyan Wang, Fan Zhang, Jialan Liang, Qian Wang, Zuo Luan, Yinxiang Yang

Institutions: Beijing Children’s Hospital (China); Chinese PLA General Hospital (China); Chinese People's Liberation Army (China)

Journal: Cell death discovery, volume 12, issue 1, article 112

Dates: received 12 June 2025; accepted 16 February 2026; published online 2 March 2026

Type: Review · Language: English

License: CC BY

Identifiers: DOI 10.1038/s41420-026-02971-w · PMID 41771834 · PMCID PMC12979777 · OpenAlex W7133221785

Open access: gold, a free copy (OpenAlex)

Status: code verified

Categories: genetics / omics (modality), human (organism), cellular / molecular (subfield)

Methods: Spectral & time-frequency, Statistics, Smoothing, state filtering, decompositions

Keywords: Stem-cell differentiation, Development

Topic: Neurogenesis and neuroplasticity mechanisms (Developmental Neuroscience, Neuroscience), according to OpenAlex

Funding: National Key Researsh and Development Program of China,2017YFA0104200

Citations: cited by 2 papers (Europe PMC); 52 references in the paper

Abstract

Heterogeneity is widely recognised across different cell types. Human oligodendrocyte progenitor cells (hOPCs), essential for myelination, exhibit considerable heterogeneity, which has not been fully characterised. In the current study, by examining the transcriptome of hOPCs at the single-cell level, three distinct subclusters were identified: PRE-OPCs, OPCs, and PRE-OLs. Single-cell RNA-sequencing and RNA-Scope detected high platelet-derived growth factor receptor alpha (PDGFRA) expression. PDGFR-α+ hOPCs exhibited greater myelination, migration, and proliferation capabilities compared to both unsorted hOPCs and PDGFR-α– hOPCs. These enhanced functions may be associated with the activation of the PI3K-AKT-mTOR and TGF-β signalling pathways, which support oligodendrocyte differentiation.

hOPCs were induced by hNSCs, their characteristics were identified. RNA-Scope and single-cell RNA Seq sequencing showed PDGFRA were highly expressed at mRNA and protein level. hOPCs were sorted by MACS using PDGFR-α beads. The myelination, migration, and proliferation abilities of PDGFR-α+ hOPCs were higher than that of un-sorting hOPCs and PDGFR-α– hOPCs, possibly being associated with the activation of PI3K–AKT–mTOR and TGF-β signalling pathways, which support oligodendrocyte differentiation (Partly created with Scientific Image and Illustration Software BioRender).

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 5 matches between paragraphs and lines of code.

satijalab/seurat

License: other

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: 586015abde10618ecb32d3fe632267a83317a08d, 21 September 2026

Languages: R (114), C++ (8), C/C++ (4), C (1)

Size: 455 files, 127 scripts

Software Heritage: archived

Found in: the text, “Single-cell RNA sequencing (scRNA‑seq) of hOPCs”

Holds: README, license file, environment (DESCRIPTION), tests, continuous integration, documentation, 70 notebooks

Not found: CITATION.cff

Tools: Seurat (75 files), ggplot2 (48 files), patchwork (27 files), tidyverse (19 files), cowplot (10 files), reshape2 (3 files), SingleCellExperiment (3 files), Plotly (2 files), data.table (1 file), DESeq2 (1 file), Harmony (1 file), igraph (1 file), limma (1 file), Monocle 3 (1 file), reticulate (1 file), UMAP (1 file)

Availability: 1 check, the latest on 27 September 2026: the link answers

27 September 2026: the link answers

130 files

R/RcppExports.R, R, 123 lines
R/clustering.R, R, 1,908 lines
R/convenience.R, R, 512 lines
R/data.R, R, 55 lines
R/differential_expression.R, R, 2,569 lines, 1 match
R/dimensional_reduction.R, R, 2,935 lines
R/generics.R, R, 842 lines, 1 match
R/integration.R, R, 5,665 lines
R/integration5.R, R, 760 lines
R/mixscape.R, R, 1,347 lines
R/objects.R, R, 3,369 lines
R/preprocessing.R, R, 5,956 lines, 1 match
R/preprocessing5.R, R, 1,891 lines
R/reexports.R, R, 454 lines
R/roxygen.R, R, 61 lines
R/sketching.R, R, 803 lines
R/tree.R, R, 392 lines
R/utilities.R, R, 3,220 lines
R/visualization.R, R, 6,161 lines, 1 match
R/zzz.R, R, 110 lines
man/roxygen/templates/note-reqdpkg.R, R, 3 lines
man/roxygen/templates/param-dotsi.R, R, 1 line
man/roxygen/templates/param-dotsm.R, R, 1 line
man/roxygen/templates/section-future.R, R, 13 lines
man/roxygen/templates/section-progressr.R, R, 8 lines
man/roxygen/templates/seealso-methods.R, R, 1 line
src/ModularityOptimizer.cpp, C++, 1,012 lines
src/ModularityOptimizer.h, C/C++, 159 lines
src/RModularityOptimizer.cpp, C++, 173 lines
src/RcppExports.cpp, C++, 445 lines
src/data_manipulation.cpp, C++, 464 lines
src/data_manipulation.h, C/C++, 51 lines
src/fast_NN_dist.cpp, C++, 70 lines
src/integration.cpp, C++, 177 lines
src/integration.h, C/C++, 37 lines
src/snn.cpp, C++, 126 lines
src/snn.h, C/C++, 21 lines
src/stats.cpp, C++, 42 lines
src/valid_pointer.c, C, 6 lines
tests/testthat.R, R, 14 lines
tests/testthat/test_data_manipulation.R, R, 209 lines
tests/testthat/test_differential_expression.R, R, 704 lines
tests/testthat/test_dimensional_reduction.R, R, 190 lines
tests/testthat/test_find_clusters.R, R, 129 lines
tests/testthat/test_integratedata.R, R, 117 lines
tests/testthat/test_integration.R, R, 487 lines
tests/testthat/test_integration5.R, R, 544 lines
tests/testthat/test_load_10X.R, R, 281 lines
tests/testthat/test_modularity_optimizer.R, R, 144 lines
tests/testthat/test_objects.R, R, 25 lines
tests/testthat/test_preprocessing.R, R, 631 lines
tests/testthat/test_sketching.R, R, 291 lines
tests/testthat/test_spatial.R, R, 509 lines
tests/testthat/test_transferdata.R, R, 143 lines
tests/testthat/test_tree.R, R, 54 lines
tests/testthat/test_utilities.R, R, 476 lines
tests/testthat/test_visualization.R, R, 54 lines
vignettes/COVID_SCTMapping.Rmd, R, 180 lines
vignettes/ParseBio_sketch_integration.Rmd, R, 180 lines
vignettes/announcements.Rmd, R, 56 lines
vignettes/archive.Rmd, R, 41 lines
vignettes/atacseq_integration_vignette.Rmd, R, 231 lines
vignettes/cell_cycle_vignette.Rmd, R, 139 lines
vignettes/conversion_vignette.Rmd, R, 132 lines
vignettes/de_vignette.Rmd, R, 212 lines
vignettes/dim_reduction_vignette.Rmd, R, 119 lines
vignettes/essential_commands.Rmd, R, 427 lines
vignettes/extensions.Rmd, R, 38 lines
vignettes/future_vignette.Rmd, R, 133 lines
vignettes/get_started.Rmd, R, 159 lines
vignettes/get_started_v5.Rmd, R, 140 lines
vignettes/get_started_v5_new.Rmd, R, 172 lines
vignettes/hashing_vignette.Rmd, R, 285 lines
vignettes/install.Rmd, R, 120 lines
vignettes/install_v5.Rmd, R, 84 lines
vignettes/integration_introduction.Rmd, R, 246 lines
vignettes/integration_large_datasets.Rmd, R, 119 lines
vignettes/integration_mapping.Rmd, R, 178 lines
vignettes/integration_rpca.Rmd, R, 187 lines
vignettes/interaction_vignette.Rmd, R, 149 lines
vignettes/merge_vignette.Rmd, R, 109 lines
vignettes/mixscape_vignette.Rmd, R, 375 lines
vignettes/multimodal_reference_mapping.Rmd, R, 387 lines
vignettes/multimodal_vignette.Rmd, R, 243 lines
vignettes/pbmc3k_tutorial.Rmd, R, 376 lines
vignettes/sctransform_v2_vignette.Rmd, R, 230 lines
vignettes/sctransform_vignette.Rmd, R, 144 lines
vignettes/seurat5_archive.Rmd, R, 41 lines
vignettes/seurat5_atacseq_integration_vignette.Rmd, R, 221 lines
vignettes/seurat5_bpcells_interaction_vignette.Rmd, R, 212 lines
vignettes/seurat5_cell_cycle_vignette.Rmd, R, 150 lines
vignettes/seurat5_conversion_vignette.Rmd, R, 132 lines
vignettes/seurat5_de_vignette.Rmd, R, 143 lines
vignettes/seurat5_dim_reduction_vignette.Rmd, R, 119 lines
vignettes/seurat5_essential_commands.Rmd, R, 220 lines
vignettes/seurat5_extensions.Rmd, R, 35 lines
vignettes/seurat5_future_vignette.Rmd, R, 134 lines
vignettes/seurat5_get_started.Rmd, R, 159 lines
vignettes/seurat5_hashing_vignette.Rmd, R, 291 lines
vignettes/seurat5_integration.Rmd, R, 220 lines
vignettes/seurat5_integration_bridge.Rmd, R, 284 lines
vignettes/seurat5_integration_introduction.Rmd, R, 273 lines
vignettes/seurat5_integration_large_datasets.Rmd, R, 125 lines
vignettes/seurat5_integration_mapping.Rmd, R, 211 lines
vignettes/seurat5_integration_rpca.Rmd, R, 191 lines
vignettes/seurat5_interaction_vignette.Rmd, R, 146 lines
vignettes/seurat5_merge_vignette.Rmd, R, 109 lines
vignettes/seurat5_mixscape_vignette.Rmd, R, 372 lines
vignettes/seurat5_multimodal_reference_mapping.Rmd, R, 396 lines
vignettes/seurat5_multimodal_vignette.Rmd, R, 240 lines
vignettes/seurat5_pbmc3k_tutorial.Rmd, R, 395 lines
vignettes/seurat5_sctransform_integration.Rmd, R, 90 lines
vignettes/seurat5_sctransform_v2_vignette.Rmd, R, 231 lines
vignettes/seurat5_sctransform_vignette.Rmd, R, 156 lines
vignettes/seurat5_sketch_analysis.Rmd, R, 209 lines
vignettes/seurat5_spatial_vignette.Rmd, R, 524 lines
vignettes/seurat5_spatial_vignette_2.Rmd, R, 570 lines
vignettes/seurat5_v4_changes.Rmd, R, 38 lines
vignettes/seurat5_visualization_vignette.Rmd, R, 249 lines
vignettes/seurat5_weighted_nearest_neighbor_analysis.Rmd, R, 451 lines
vignettes/spatial_vignette.Rmd, R, 603 lines
vignettes/spatial_vignette_2.Rmd, R, 734 lines
vignettes/v4_changes.Rmd, R, 38 lines
vignettes/visiumhd_analysis_cell_segmentations.Rmd, R, 451 lines
vignettes/visiumhd_analysis_vignette.Rmd, R, 559 lines
vignettes/visualization_vignette.Rmd, R, 241 lines
vignettes/weighted_nearest_neighbor_analysis.Rmd, R, 442 lines
LICENSE, License, 2 lines
LICENSE.md, License, 21 lines
README.md, Text, 23 lines

At the source: github.com/satijalab/seurat

chris-mcginnis-ucsf/DoubletFinder

License: none: the authors keep all their rights

State: the link answers, verified on 29 September 2026

Evidence: files inventoried

Commit: 1b244d8f0d54b4b1cb4365639931bbb16f01e1cd, 21 March 2025

Languages: R (12)

Size: 38 files, 12 scripts

Software Heritage: not archived

Found in: the text, “Single-cell RNA sequencing (scRNA‑seq) of hOPCs”

Holds: README, environment (DESCRIPTION), tests, documentation

Not found: license file, CITATION.cff, continuous integration

Tools: Seurat (4 files)

Availability: 1 check, the latest on 29 September 2026: the link answers

29 September 2026: the link answers

13 files

R/bimodality_coefficient.R, R, 19 lines
R/data.R, R, 12 lines
R/doubletFinder.R, R, 209 lines, 1 match
R/find.pK.R, R, 89 lines
R/kurtosis.R, R, 18 lines
R/modelHomotypic.R, R, 41 lines
R/parallel_paramSweep.R, R, 138 lines
R/paramSweep.R, R, 110 lines
R/skewness.R, R, 17 lines
R/summarizeSweep.R, R, 85 lines
tests/testthat.R, R, 5 lines
tests/testthat/test-doubletFinder_v3.R, R, 29 lines
README.md, Text, 202 lines

At the source: github.com/chris-mcginnis-ucsf/DoubletFinder

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
139 scripts, each with its path and the digest of its content;
5 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability

Raw data from scRNA-seq have been deposited at GSA-human. Raw data reported in this paper will be shared by the lead contact upon request. Raw data of scRNA-seq at [https://ngdc.cncb.ac.cn/gsa-human/s/FdDsn5B5] and [https://ngdc.cncb.ac.cn/gsa-human/s/z66g0yYU].

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 30 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 11 authors, 2 keywords, 1 funder, 51 references.

Cite

This paper

Ye, D., Zhou, H., Qu, S., Wang, Z., Zhang, F., Wang, X., Zhao, Y., Liang, J., Wang, Q., Luan, Z., & Yang, Y. (2026). Multi-omics reveals heterogeneity and functional populations of oligodendrocyte progenitor cells induced by human neural stem cells. Cell death discovery, 12(1), 112. https://doi.org/10.1038/s41420-026-02971-w

BibTeX

@article{ye2026multi, author = {Ye, Dou and Zhou, Haipeng and Qu, Suqing and Wang, Zhaoyan and Zhang, Fan and Wang, Xiaohua and Zhao, Yuan and Liang, Jialan and Wang, Qian and Luan, Zuo and Yang, Yinxiang}, title = {{Multi-omics reveals heterogeneity and functional populations of oligodendrocyte progenitor cells induced by human neural stem cells}}, journal = {Cell death discovery}, year = {2026}, month = mar, volume = {12}, number = {1}, pages = {112}, publisher = {Nature Publishing Group}, issn = {2058-7716}, doi = {10.1038/s41420-026-02971-w}, url = {https://doi.org/10.1038/s41420-026-02971-w}, pmid = {41771834}, pmcid = {PMC12979777} }

RIS

TY - JOUR AU - Ye, Dou AU - Zhou, Haipeng AU - Qu, Suqing AU - Wang, Zhaoyan AU - Zhang, Fan AU - Wang, Xiaohua AU - Zhao, Yuan AU - Liang, Jialan AU - Wang, Qian AU - Luan, Zuo AU - Yang, Yinxiang TI - Multi-omics reveals heterogeneity and functional populations of oligodendrocyte progenitor cells induced by human neural stem cells T2 - Cell death discovery J2 - Cell Death Discov PY - 2026 DA - 2026/03/02 VL - 12 IS - 1 SP - 112 SN - 2058-7716 PB - Nature Publishing Group DO - 10.1038/s41420-026-02971-w UR - https://doi.org/10.1038/s41420-026-02971-w LA - en ER -

CSL-JSON

{ "id": "10.1038/s41420-026-02971-w", "type": "article-journal", "title": "Multi-omics reveals heterogeneity and functional populations of oligodendrocyte progenitor cells induced by human neural stem cells", "container-title": "Cell death discovery", "author": [ { "family": "Ye", "given": "Dou" }, { "family": "Zhou", "given": "Haipeng" }, { "family": "Qu", "given": "Suqing" }, { "family": "Wang", "given": "Zhaoyan" }, { "family": "Zhang", "given": "Fan" }, { "family": "Wang", "given": "Xiaohua" }, { "family": "Zhao", "given": "Yuan" }, { "family": "Liang", "given": "Jialan" }, { "family": "Wang", "given": "Qian" }, { "family": "Luan", "given": "Zuo" }, { "family": "Yang", "given": "Yinxiang" } ], "container-title-short": "Cell Death Discov", "volume": "12", "issue": "1", "page": "112", "DOI": "10.1038/s41420-026-02971-w", "PMID": "41771834", "PMCID": "PMC12979777", "ISSN": "2058-7716", "publisher": "Nature Publishing Group", "URL": "https://doi.org/10.1038/s41420-026-02971-w", "language": "en", "issued": { "date-parts": [ [ 2026, 3, 2 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.7554/elife.93640 [code]: Sibling chimerism among microglia in marmosets.
Journal: eLife
In common: Monocle 3, Harmony, SingleCellExperiment, 13 other tools, genetics / omics, cellular / molecular
[2] doi:10.1186/s12974-026-03838-8 [code]: Acarbose modulates microglial Pkm2 acetylation to reshape immunometabolism and preserve retinal neurons after ischemia-reperfusion.
Journal: Journal of neuroinflammation
In common: Monocle 3, Harmony, SingleCellExperiment, 13 other tools, genetics / omics, cellular / molecular
[3] doi:10.1002/ctm2.70683 [code]: Niacin promotes motor function recovery after spinal cord injury via Hcar2-dependent microglia immunometabolic regulation.
Journal: Clinical and translational medicine
In common: Monocle 3, Harmony, SingleCellExperiment, 13 other tools, genetics / omics, cellular / molecular
[4] doi:10.1038/s44318-026-00818-9 [code]: FAM134B-mediated ER-phagy degrades APP and suppresses Alzheimer's disease pathology.
Journal: The EMBO journal
In common: Monocle 3, Harmony, SingleCellExperiment, 13 other tools, cellular / molecular
[5] doi:10.1016/j.xcrm.2026.102766 [code]: A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: Monocle 3, Harmony, SingleCellExperiment, 12 other tools, genetics / omics, cellular / molecular, 1 reference
[6] doi:10.1038/s41467-026-71803-3 [code]: Charting the transition from in vitro gliogenesis to the in vivo maturation of human glial progenitor cells transplanted into the hypomyelinated mouse brain.
Journal: Nature communications
In common: reticulate, DESeq2, Plotly, 6 other tools, genetics / omics, cellular / molecular, 7 references
[7] doi:10.1038/s41593-026-02367-0 [code]: A reproducible three-dimensional model of human brain tissue to investigate physiological and disease-associated microglia phenotypes.
Journal: Nature neuroscience
In common: Monocle 3, Harmony, SingleCellExperiment, 11 other tools, cellular / molecular, 1 reference
[8] doi:10.3390/ijms27104466 [code]: Uncovering the Key Circuit FOSL2/FOS/EGR3/EGR1, Contributing to the Hyperexcitability of Excitatory Neurons in the Epileptic Temporal Cortex and Hippocampus.
Journal: International journal of molecular sciences
In common: Monocle 3, reticulate, limma, 10 other tools, genetics / omics, 2 references
[9] doi:10.1002/imt2.70163 [code]: Spatial multi-omics unveils sphingolipid metabolic reprogramming within the retinal pathological niche.
Journal: iMeta
In common: Monocle 3, Harmony, SingleCellExperiment, 11 other tools, genetics / omics, cellular / molecular
[10] doi:10.1016/j.xcrm.2026.102651 [code]: Integrative CSF profiling identifies disease-specific immune responses in leptomeningeal disease.
Journal: Cell reports. Medicine
In common: Harmony, SingleCellExperiment, reticulate, 11 other tools, genetics / omics, cellular / molecular

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.