#' @importFrom ggplot2 ggproto
NULL
####
# Spatial plots ####
####
####
## Spatial Identity Plot ####
####
#' vrSpatialPlot
#'
#' Plotting identification of spatially resolved cells, spots, and ROI on associated images from multiple assays in a VoltRon object.
#'
#' @param object a VoltRon object
#' @param group.by a column of metadata from \link{Metadata} used as grouping label for the spatial entities
#' @param plot.segments plot segments from \link{vrSegments} instead of points
#' @param group.ids a subset of categories defined in metadata column from \code{group.by}
#' @param colors the color set for group.by. Should be of the same size of group.id (if specified) or unique elements in group.by
#' @param n.tile should points be aggregated into tiles before visualization (see \link{geom_tile}). Applicable only for cells and molecules
#' @param assay assay name (exp: Assay1) or assay class (exp: Visium, Xenium), see \link{SampleMetadata}.
#' if NULL, the default assay will be used, see \link{vrMainAssay}.
#' @param graph.name if not NULL, the spatial graph is with name \code{graph.name} is visualized as well, see \link{vrGraphNames}
#' @param graph.edge.color the colors of the graph edges, if \code{graph.name} is not NULL.
#' @param reduction used by \code{vrSpatialPlotVitessce} to visualize an embedding alongside with the spatial plot.
#' @param ncol column wise number of plots, for \link{ggarrange}
#' @param nrow row wise number of plots, for \link{ggarrange}
#' @param font.size font size
#' @param pt.size point size
#' @param cell.shape the shape of the points representing cells, see \link{geom_point}
#' @param alpha alpha level of colors of visualized points and segments
#' @param label if TRUE, the labels of the ROI assays will be visualized
#' @param spatial the name of the main spatial system
#' @param channel the name of the channel associated with the image
#' @param background.color the color of plot background if a channel is not specified, or the spatial coord system doesnt have an image.
#' @param background (DEPRECATED) the background of the plot. Either an image name, see \link{vrImageNames} or a vector of length two with image name
#' and a channel name, see \link{vrImageChannelNames}. Type "black" or "white" for black or white backgrounds. if NULL, the main image (\link{vrMainSpatial})
#' and main channel (\link{vrMainChannel}) will be in the background. Otherwise the background will be grey.
#' @param reg TRUE if registered coordinates of the main image (\link{vrMainSpatial}) is requested
#' @param crop whether to crop an image of a spot assay to the extend of spots
#' @param legend.pt.size the size of points at the legend
#' @param legend.text.size the size of the text at the legend
#' @param scale.image if TRUE, background image will be scaled down to a low resolution (width: 1000px)
#' @param legend.loc the location of the legend, default is "right"
#' @param common.legend whether to use a common legend for all plots, see \link{ggarrange}
#' @param collapse.plots whether to combine all ggplots
#' @param interactive if TRUE, run interactive plot
#' @param shiny.options a list of shiny options (browser, host, port etc.) passed \code{options} arguement of \link{shinyApp}. For more information, see \link{runApp}
#'
#' @import ggplot2
#' @importFrom ggpubr ggarrange
#'
#' @export
vrSpatialPlot <- function(object, group.by = "Sample", plot.segments = FALSE, group.ids = NULL, colors = NULL, n.tile = 0,
assay = NULL, graph.name = NULL, graph.edge.color = "orange", reduction = NULL, ncol = 2, nrow = NULL, font.size = 2, pt.size = 2,
cell.shape = 21, alpha = 1, label = FALSE, spatial = NULL, channel = NULL, background.color = NULL,
background = NULL, reg = FALSE, crop = FALSE, legend.pt.size = 2, legend.text.size = 14,
scale.image = TRUE, legend.loc = "right", common.legend = TRUE, collapse.plots = TRUE, interactive = FALSE,
shiny.options = list()) {
# check object for zarr
if(is.character(object)){
if(grepl(".zarr$", object)){
return(vrSpatialPlotVitessce(zarr.file = object, group.by = group.by, reduction = reduction,
shiny.options = shiny.options))
}
}
# check object
if(!inherits(object, "VoltRon"))
stop("Please provide a VoltRon object!")
# sample metadata
sample.metadata <- SampleMetadata(object)
# get assay names
assay_names <- vrAssayNames(object, assay = assay)
# interactive plotting
if(interactive){
if(length(assay_names) > 1){
stop("Only one assay can be visualized with the interactive plot at a time.")
} else{
gg <- vrSpatialPlot(object, group.by = group.by, plot.segments = plot.segments, group.ids = group.ids, colors = colors, n.tile = n.tile,
assay = assay, graph.name = graph.name, reduction = reduction, ncol = ncol, nrow = nrow, font.size = font.size,
pt.size = pt.size, cell.shape = cell.shape, alpha = alpha, label = label, spatial = spatial, channel = channel,
background.color = background.color, background = background, reg = reg, crop = crop, legend.pt.size = legend.pt.size,
legend.text.size = legend.text.size, scale.image = FALSE, legend.loc = legend.loc, common.legend = common.legend,
collapse.plots = collapse.plots, interactive = FALSE)
return(vrSpatialPlotInteractive(plot_g = gg, shiny.options = shiny.options))
}
}
# get entity type and metadata
metadata <- Metadata(object, assay = assay)
# configure titles
plot_title <- as.list(apply(sample.metadata[assay_names,], 1, function(x) paste(x["Sample"], x["Layer"], sep = ", ")))
names(plot_title) <- assay_names
# adjust group.ids
if(!group.by %in% colnames(metadata))
stop("The column '", group.by, "' was not found in the metadata!")
if(is.null(group.ids)){
group.ids <- unique(metadata[[group.by]])
}
# check group.id
levels_group.ids <- as.character(group.ids)
if(all(!is.na(suppressWarnings(as.numeric(levels_group.ids))))){
levels_group.ids <- sort(as.numeric(levels_group.ids))
}
group.ids <- factor(group.ids, levels = levels_group.ids)
# check colors
if(!is.null(colors)){
if(length(colors) != length(levels_group.ids)){
stop("Please provide colors whose length is equal to group.ids (if specified) or unique elements of group.by")
}
if(is.null(names(colors))){
stop("Name of each color has to be specified given as names(colors)")
}
} else{
colors <- hue_pal(length(group.ids))
names(colors) <- group.ids
}
# for each assay
i <- 1
gg <- list()
for(assy in assay_names){
# get assay
cur_assay <- object[[assy]]
cur_metadata <- subset_metadata(metadata, assays = assy)
# get graph
if(!is.null(graph.name)){
if(graph.name %in% vrGraphNames(object)){
graph <- vrGraph(object, assay = assy, graph.type = graph.name)
} else {
stop("the graph with named '", graph.name, "' was not found in the graph list!")
}
} else {
graph <- NULL
}
# check group.by
levels_group.by <- as.character(unique(metadata[[group.by]][!is.na(metadata[[group.by]])]))
if(all(!is.na(suppressWarnings(as.numeric(levels_group.by))))){
levels_group.by <- sort(as.numeric(levels_group.by))
}
cur_metadata[[group.by]] <- factor(cur_metadata[[group.by]], levels = levels_group.by)
# visualize
p_title <- plot_title[[assy]]
gg[[i]] <- vrSpatialPlotSingle(assay = cur_assay, metadata = cur_metadata,
group.by = group.by, plot.segments = plot.segments, group.ids = group.ids, colors = colors, n.tile = n.tile, graph = graph, graph.edge.color = graph.edge.color, font.size = font.size,
pt.size = pt.size, alpha = alpha, cell.shape = cell.shape, plot_title = p_title, background = background, spatial = spatial, channel = channel, background.color = background.color, reg = reg,
crop = crop, legend.pt.size = legend.pt.size, legend.text.size = legend.text.size, scale.image = scale.image)
i <- i + 1
}
# return a list of plots or a single one
if(collapse.plots){
if(length(assay_names) > 1){
if(length(gg) < ncol) ncol <- length(gg)
return(ggpubr::ggarrange(plotlist = gg, ncol = ncol, nrow = ceiling(length(gg)/ncol), common.legend = common.legend, legend = legend.loc))
} else {
return(gg[[1]] + theme(legend.position=legend.loc))
}
} else {
for(i in seq_len(length(assay_names))){
gg[[i]] <- gg[[i]] + theme(legend.position=legend.loc)
}
return(gg)
}
}
#' vrSpatialPlotSingle
#'
#' Plotting a single assay from a VoltRon object. We plot the identification of spatially resolved cells, spots, and ROI on associated images.
#'
#' @param assay vrAssay object
#' @param metadata the metadata associated with the assay
#' @param group.by a column of metadata from \link{Metadata} used as grouping label for the spatial entities
#' @param plot.segments plot segments from \link{vrSegments} instead of points
#' @param group.ids a subset of categories defined in metadata column from \code{group.by}
#' @param colors the color set for group.by. Should be of the same size of group.id (if specified) or unique elements in group.by
#' @param n.tile should points be aggregated into tiles before visualization (see \link{geom_tile}). Applicable only for cells and molecules
#' @param graph if not NULL, the graph is added to the plot
#' @param graph.edge.color the color of graph edges, if \code{graph} is not NULL.
#' @param font.size font sizes
#' @param pt.size point size
#' @param cell.shape the shape of the points representing cells, see \link{geom_point}
#' @param alpha alpha level of colors of visualized points and segments
#' @param plot_title the title of the single plot
#' @param spatial the name of the main spatial system
#' @param channel the name of the channel associated with the image
#' @param background.color the color of plot background if a channel is not specified, or the spatial coord system doesnt have an image.
#' @param background (DEPRECATED) the background of the plot. Either an image name, see \link{vrImageNames} or a vector of length two with image name
#' and a channel name, see \link{vrImageChannelNames}. Type "black" or "white" for black or white backgrounds. if NULL, the main image (\link{vrMainSpatial})
#' and main channel (\link{vrMainChannel}) will be in the background. Otherwise the background will be grey.
#' @param reg TRUE if registered coordinates of the main image (\link{vrMainSpatial}) is requested
#' @param crop whether to crop an image of a spot assay to the extend of spots
#' @param legend.pt.size the size of points at the legend
#' @param legend.text.size the size of the text at the legend
#' @param scale.image if TRUE, background image will be scaled down to a low resolution (width: 1000px)
#'
#' @import ggplot2
#' @importFrom igraph get.data.frame
#' @importFrom magick image_info image_read
#'
#' @noRd
vrSpatialPlotSingle <- function(assay, metadata, group.by = "Sample", plot.segments = FALSE, group.ids = NULL, colors = NULL, n.tile = 0,
graph = NULL, graph.edge.color = "orange", font.size = 2, pt.size = 2, cell.shape = 16, alpha = 1, plot_title = NULL, spatial = NULL,
channel = NULL, background.color = NULL, background = NULL, reg = FALSE, crop = FALSE, legend.pt.size = 2,
legend.text.size = 14, scale.image = TRUE){
# plot
g <- ggplot()
# add assay to ggplot
g$voltron_params <- list()
g$voltron_params$assay <- vrAssayNames(assay)
# add image and background
image <- vrSpatialPlotImage(g, assay, background, scale.image, spatial = spatial,
channel = channel, background.color = background.color)
g <- image$plot
info <- image$info
background.color <- image$background.color
g$voltron_params$spatial_name <- spatial_name <- image$spatial
g$voltron_params$scale_factors <- scale_factors <- image$scale_factors
# coords
coords <- vrCoordinates(assay, spatial_name = spatial_name, reg = reg)
if(!inherits(coords, "IterableMatrix")){
coords <- as.data.frame(coords)
}
coords <- coords/scale_factors
segments <- vrSegments(assay, spatial_name = spatial_name)
# plotting features
if(!group.by %in% colnames(metadata))
stop("The column '", group.by, "' was not found in the metadata!")
if(inherits(metadata, "data.table")){
cur_group.by <- metadata[,get(names(metadata)[which(colnames(metadata) == group.by)])]
names(cur_group.by) <- metadata$id
} else {
cur_group.by <- metadata[,group.by]
if(!is.null(rownames(metadata))){
names(cur_group.by) <- rownames(metadata)
} else {
names(cur_group.by) <- as.vector(metadata$id)
}
}
# group ids
if(!is.null(group.ids)){
len_set_diff <- length(setdiff(group.ids, cur_group.by))
if(len_set_diff > 0){
# warning("Some groups defined in group.ids does not exist in group.by!")
} else if(len_set_diff == length(group.ids)){
stop("None of the groups defined in group.ids exist in group.by!")
} else {
segments <- segments[cur_group.by %in% group.ids]
}
cur_group.by <- droplevels(cur_group.by[cur_group.by %in% group.ids])
coords <- coords[names(cur_group.by),]
}
# change levels of groups
cur_group.by <- factor(cur_group.by, levels = group.ids)
# add groups
coords <- as.matrix(coords)
coords <- data.frame(coords, cur_group.by[rownames(coords), drop = FALSE])
colnames(coords)[length(colnames(coords))] <- group.by
# set up the limits
g <- vrSpatialExtent(g, assay, coords, crop, info)
# set colors
names_colors <- factor(names(colors))
# visualize based on points type
if(vrAssayTypes(assay) == "ROI"){
polygon_data <- NULL
circle_data <- NULL
for(i in seq_len(length(segments))){
if(nrow(segments[[i]]) > 1){
cur_segment <- segments[[i]][,c("x","y")]
suppressWarnings({
cur_data <- as.data.frame(cbind(cur_segment, names(segments)[i], cur_group.by[i]))
})
colnames(cur_data) <- c("x", "y", "segment", "group.by")
cur_data[,c("x", "y")] <- cur_data[,c("x", "y")]/scale_factors
polygon_data <- as.data.frame(rbind(polygon_data, cur_data))
} else {
cur_segment <- segments[[i]][,c("x","y", "rx", "ry")]
suppressWarnings({
cur_data <- as.data.frame(cbind(cur_segment, names(segments)[i], cur_group.by[i]))
})
colnames(cur_data) <- c("x", "y", "rx", "ry", "segment", "group.by")
cur_data[,c("x", "y","rx", "ry")] <- cur_data[,c("x", "y","rx", "ry")]/scale_factors
circle_data <- as.data.frame(rbind(circle_data, cur_data))
}
}
if(!is.null(polygon_data)){
g <- g +
geom_polygon(aes(x = x, y = y, fill = group.by, group = segment), data = polygon_data, alpha = alpha, show.legend = TRUE)
}
if(!is.null(circle_data)){
if(!requireNamespace('ggforce'))
stop("Please install ggforce package!: install.packages('ggforce')")
g <- g +
ggforce::geom_ellipse(aes(x0 = as.numeric(x), y0 = as.numeric(y), a = as.numeric(rx), b = as.numeric(ry), angle = 0,
fill = group.by, group = segment), data = circle_data, lwd = 0, alpha = alpha, show.legend = TRUE)
}
g <- g +
scale_fill_manual(values = colors, labels = names_colors, drop = FALSE, limits = names_colors,
guide = guide_legend(override.aes = list(alpha = 1))) +
guides(fill = guide_legend(title = group.by))
# spot visualization
} else if(vrAssayTypes(assay) == "spot"){
# rasterize if requested or needed
if(n.tile > 0 || nrow(coords) > 50000){
if(n.tile == 0)
n.tile <- 1000
g <- vrGroupPlotTiling(g = g, data = coords, group.by = group.by, n.tile = n.tile, alpha = alpha)
} else {
spot.type <- vrAssayParams(assay, param = "spot.type")
spot.type <- ifelse(is.null(spot.type), "circle", spot.type)
g <- g +
geom_spot(mapping = aes(x = .data[["x"]], y = .data[["y"]], fill = .data[[group.by]]), coords, shape = 21, alpha = alpha, show.legend = TRUE,
spot.radius = vrAssayParams(assay, param = "vis.spot.radius")/scale_factors,
spot.type = spot.type)
# add if a graph exists
if(!is.null(graph)){
g <- g + addGraph(graph = graph, coords = coords, background = graph.edge.color)
}
}
# style, color and text
g <- g +
scale_fill_manual(values = colors, labels = names_colors, drop = FALSE, limits = names_colors) +
guides(fill = guide_legend(override.aes=list(shape = 21, size = 4, lwd = 0.1)))
# cell visualization
} else if(vrAssayTypes(assay) %in% c("cell", "tile")) {
if(plot.segments){
if(length(segments) == 0) {
stop("No Segments are available in this assay!")
} else {
polygon_data <- do.call(rbind,segments)
polygon_data[,c("x", "y")] <- polygon_data[,c("x", "y")]/scale_factors
len_segments <- vapply(segments, nrow, numeric(1))
polygon_data <- data.frame(polygon_data, segment = rep(names(segments), len_segments), group.by = rep(cur_group.by, len_segments))
g <- g +
geom_polygon(aes(x = .data[["x"]], y = .data[["y"]], fill = .data[["group.by"]], group = segment), data = polygon_data, alpha = alpha, show.legend = TRUE) +
scale_fill_manual(values = colors, labels = names_colors, drop = FALSE, limits = names_colors, name = group.by, guide = guide_legend(order = 1))
}
} else {
# rasterize if requested or needed
if(n.tile > 0 || nrow(coords) > 50000){
if(n.tile == 0)
n.tile <- 1000
g <- vrGroupPlotTiling(g = g, data = coords, group.by = group.by, n.tile = n.tile, alpha = alpha)
} else {
g <- g +
geom_point(mapping = aes(x = .data[["x"]], y = .data[["y"]], fill = .data[[group.by]], color = .data[[group.by]]),
coords, shape = cell.shape, size = rel(pt.size), alpha = alpha, show.legend = TRUE)
}
# style, color and text
g <- g +
scale_fill_manual(values = colors, labels = names_colors, drop = FALSE, limits = names_colors, name = group.by, guide = guide_legend(order = 1)) +
scale_color_manual(values = colors, labels = names_colors, drop = FALSE, limits = names_colors, name = group.by, guide = guide_legend(order = 1)) +
theme(legend.text=element_text(size=legend.text.size), legend.title=element_text(size=legend.text.size))
# add if a graph exists
if(!is.null(graph)){
g <- g + addGraph(graph = graph, coords = coords, background = graph.edge.color)
}
}
} else if(vrAssayTypes(assay) == "molecule") {
# rasterize if requested or needed
if(n.tile > 0 || nrow(coords) > 50000){
if(n.tile == 0)
n.tile <- 1000
g <- vrGroupPlotTiling(g = g, data = coords, group.by = group.by, n.tile = n.tile, alpha = alpha)
} else {
g <- g +
geom_point(mapping = aes(x = .data[["x"]], y = .data[["y"]], fill = .data[[group.by]], color = .data[[group.by]]),
coords, shape = cell.shape, size = rel(pt.size), alpha = alpha, show.legend = TRUE)
}
g <- g +
scale_fill_manual(values = colors, labels = names_colors, drop = FALSE, limits = names_colors, name = group.by, guide = guide_legend(order = 1)) +
scale_color_manual(values = colors, labels = names_colors, drop = FALSE, limits = names_colors, name = group.by, guide = guide_legend(order = 1)) +
theme(legend.text=element_text(size=legend.text.size), legend.title=element_text(size=legend.text.size))
} else {
stop("Only ROIs, spots, cells, molecules and tiles can be visualized with vrSpatialPlot!")
}
# more visualization parameters
g <- g +
ggtitle(plot_title) + theme(plot.title = element_text(hjust = 0.5, margin=margin(0,0,0,0), size = legend.text.size),
panel.grid.minor = element_blank(), panel.grid.major = element_blank(),
axis.line=element_blank(),axis.text.x=element_blank(), axis.text.y=element_blank(),
axis.ticks=element_blank(), axis.title.x=element_blank(), axis.title.y=element_blank(),
legend.margin = margin(0,0,0,0))
# return data
return(g)
}
#' addSpatialLayer
#'
#' adding additional layers of spatial plots to an existing \link{vrSpatialPlot}.
#'
#' @param g ggplot object
#' @param object a VoltRon object
#' @param assay assay name (exp: Assay1) or assay class (exp: Visium, Xenium), see \link{SampleMetadata}.
#' if NULL, the default assay will be used, see \link{vrMainAssay}.
#' @param group.by a column of metadata from \link{Metadata} used as grouping label for the spatial entities
#' @param plot.segments plot segments from \link{vrSegments} instead of points
#' @param group.ids a subset of categories defined in metadata column from \code{group.by}
#' @param reg TRUE if registered coordinates of the main image (\link{vrMainSpatial}) is requested
#' @param colors the color set for group.by. Should be of the same size of group.id (if specified) or unique elements in group.by
#' @param alpha alpha level of colors of visualized points and segments
#' @param n.tile should points be aggregated into tiles before visualization (see \link{geom_tile}). Applicable only for cells and molecules
#' @param pt.size point size
#' @param cell.shape the shape of the points representing cells, see \link{geom_point}
#' @param graph if not NULL, the graph is added to the plot
#' @param graph.edge.color the color of graph edges, if \code{graph} is not NULL.
#' @param spatial the name of the main spatial system
#'
#' @import ggplot2
#'
#' @export
addSpatialLayer <- function(g, object, assay, group.by = "Sample", plot.segments = FALSE, group.ids = NULL, reg = FALSE, colors = NULL, alpha = 1,
n.tile = 0, pt.size = 2, cell.shape = 21, graph = NULL, graph.edge.color = "orange", spatial = NULL){
# check package
if(!requireNamespace("ggnewscale")){
stop("Please install ggnewscale package for using multilayer visualization!: install_packages('ggnewscale')")
}
# sample metadata
sample.metadata <- SampleMetadata(object)
spatial.names <- vrSpatialNames(object, assay = "all")
# get assay names
assay_names <- vrAssayNames(object, assay = assay)
# objects and parameters
assay <- object[[assay_names]]
adj_assay <- g$voltron_params$assay
spatial_name <- if(is.null(spatial)) g$voltron_params$spatial_name else spatial
scale_factors <- g$voltron_params$scale_factors
# check adjacency
adj_sample <- unique(sample.metadata[c(assay_names, adj_assay),"Sample"])
adj_spatial <- unique(spatial.names[c(assay_names, adj_assay),"Main"])
if(length(adj_sample) != 1 || length(adj_sample) != 1){
message("Assays in addSpatialLayer() are not adjacent to vrSpatialPlot()!")
return(g)
}
# get entity type and metadata
metadata <- Metadata(object, assay = assay_names)
# coords
coords <- vrCoordinates(assay, spatial_name = spatial_name, reg = reg)
if(!inherits(coords, "IterableMatrix")){
coords <- as.data.frame(coords)
}
coords <- coords/scale_factors
segments <- vrSegments(assay, spatial_name = spatial_name)
# adjust group.ids
if(!group.by %in% colnames(metadata))
stop("The column '", group.by, "' was not found in the metadata!")
if(is.null(group.ids)){
group.ids <- unique(metadata[[group.by]])
}
# check group.id
levels_group.ids <- as.character(group.ids)
if(all(!is.na(suppressWarnings(as.numeric(levels_group.ids))))){
levels_group.ids <- sort(as.numeric(levels_group.ids))
}
group.ids <- factor(group.ids, levels = levels_group.ids)
# check colors
if(!is.null(colors)){
if(length(colors) != length(levels_group.ids)){
stop("Please provide colors whose length is equal to group.ids (if specified) or unique elements of group.by")
}
if(is.null(names(colors))){
stop("Name of each color has to be specified given as names(colors)")
}
} else{
colors <- hue_pal(length(group.ids))
names(colors) <- group.ids
}
# check group.by
levels_group.by <- as.character(unique(metadata[[group.by]][!is.na(metadata[[group.by]])]))
if(all(!is.na(suppressWarnings(as.numeric(levels_group.by))))){
levels_group.by <- sort(as.numeric(levels_group.by))
}
metadata[[group.by]] <- factor(metadata[[group.by]], levels = levels_group.by)
# plotting features
if(!group.by %in% colnames(metadata))
stop("The column '", group.by, "' was not found in the metadata!")
if(inherits(metadata, "data.table")){
cur_group.by <- metadata[,get(names(metadata)[which(colnames(metadata) == group.by)])]
names(cur_group.by) <- metadata$id
} else {
cur_group.by <- metadata[,group.by]
if(!is.null(rownames(metadata))){
names(cur_group.by) <- rownames(metadata)
} else {
names(cur_group.by) <- as.vector(metadata$id)
}
}
# group ids
if(!is.null(group.ids)){
len_set_diff <- length(setdiff(group.ids, cur_group.by))
if(len_set_diff > 0){
# warning("Some groups defined in group.ids does not exist in group.by!")
} else if(len_set_diff == length(group.ids)){
stop("None of the groups defined in group.ids exist in group.by!")
} else {
segments <- segments[cur_group.by %in% group.ids]
}
cur_group.by <- droplevels(cur_group.by[cur_group.by %in% group.ids])
coords <- coords[names(cur_group.by),]
}
# change levels of groups
cur_group.by <- factor(cur_group.by, levels = group.ids)
# add groups
coords <- as.matrix(coords)
coords <- data.frame(coords, cur_group.by[rownames(coords), drop = FALSE])
colnames(coords)[length(colnames(coords))] <- group.by
# set colors
names_colors <- factor(names(colors))
# set new scales
g <- g + ggnewscale::new_scale_color()
g <- g + ggnewscale::new_scale_fill()
# visualize based on points type
if(vrAssayTypes(assay) == "ROI"){
polygon_data <- NULL
circle_data <- NULL
for(i in seq_len(length(segments))){
if(nrow(segments[[i]]) > 1){
cur_segment <- segments[[i]][,c("x","y")]
suppressWarnings({
cur_data <- as.data.frame(cbind(cur_segment, names(segments)[i], cur_group.by[i]))
})
colnames(cur_data) <- c("x", "y", "segment", "group.by")
cur_data[,c("x", "y")] <- cur_data[,c("x", "y")]/scale_factors
polygon_data <- as.data.frame(rbind(polygon_data, cur_data))
} else {
cur_segment <- segments[[i]][,c("x","y", "rx", "ry")]
suppressWarnings({
cur_data <- as.data.frame(cbind(cur_segment, names(segments)[i], cur_group.by[i]))
})
colnames(cur_data) <- c("x", "y", "rx", "ry", "segment", "group.by")
cur_data[,c("x", "y","rx", "ry")] <- cur_data[,c("x", "y","rx", "ry")]/scale_factors
circle_data <- as.data.frame(rbind(circle_data, cur_data))
}
}
if(!is.null(polygon_data)){
g <- g +
geom_polygon(aes(x = .data$x, y = .data$y, fill = .data[["group.by"]], group = .data[["segment"]]), data = polygon_data, alpha = alpha, col = "grey")
}
if(!is.null(circle_data)){
if(!requireNamespace('ggforce'))
stop("Please install ggforce package!: install.packages('ggforce')")
g <- g +
ggforce::geom_ellipse(aes(x0 = as.numeric(x), y0 = as.numeric(y), a = as.numeric(rx), b = as.numeric(ry), angle = 0,
fill = group.by, color = "grey", group = segment), data = circle_data, lwd = 0, alpha = alpha)
}
g <- g +
scale_fill_manual(values = colors, labels = names_colors, drop = FALSE, limits = names_colors, name = group.by,
guide = guide_legend(override.aes = list(alpha = 1), order = 2)) +
scale_color_manual(values = colors, labels = names_colors, drop = FALSE, limits = names_colors, name = group.by,
guide = guide_legend(override.aes = list(alpha = 1), order = 2))
} else if(vrAssayTypes(assay) %in% c("cell", "tile")) {
if(plot.segments){
if(length(segments) == 0) {
stop("No Segments are available in this assay!")
} else {
polygon_data <- do.call(rbind,segments)
polygon_data[,c("x", "y")] <- polygon_data[,c("x", "y")]/scale_factors
len_segments <- vapply(segments, nrow, numeric(1))
polygon_data <- data.frame(polygon_data, segment = rep(names(segments), len_segments), group.by = rep(cur_group.by, len_segments))
g <- g +
geom_polygon(aes(x = .data[["x"]], y = .data[["y"]], fill = .data[["group.by"]], group = segment), data = polygon_data, alpha = alpha) +
scale_fill_manual(values = colors, labels = names_colors, drop = FALSE, limits = names_colors, name = group.by, guide = guide_legend(order = 2))
}
} else {
# rasterize if requested or needed
if(n.tile > 0 || nrow(coords) > 50000){
if(n.tile == 0)
n.tile <- 1000
g <- vrGroupPlotTiling(g = g, data = coords, group.by = group.by, n.tile = n.tile, alpha = alpha)
} else {
g <- g +
geom_point(mapping = aes(x = .data[["x"]], y = .data[["y"]], fill = .data[[group.by]], color = .data[[group.by]]),
coords, shape = cell.shape, size = rel(pt.size), alpha = alpha)
}
# style, color and text
g <- g +
scale_fill_manual(values = colors, labels = names_colors, drop = FALSE, limits = names_colors, name = group.by, guide = guide_legend(order = 2)) +
scale_color_manual(values = colors, labels = names_colors, drop = FALSE, limits = names_colors, name = group.by, guide = guide_legend(order = 2))
# add if a graph exists
if(!is.null(graph)){
g <- g + addGraph(graph = graph, coords = coords, background = graph.edge.color)
}
}
} else if(vrAssayTypes(assay) == "molecule") {
# rasterize if requested or needed
if(n.tile > 0 || nrow(coords) > 50000){
if(n.tile == 0)
n.tile <- 1000
g <- vrGroupPlotTiling(g = g, data = coords, group.by = group.by, n.tile = n.tile, alpha = alpha)
} else {
g <- g +
geom_point(mapping = aes(x = .data[["x"]], y = .data[["y"]], fill = .data[[group.by]], color = .data[[group.by]]),
coords, shape = cell.shape, size = rel(pt.size), alpha = alpha)
}
g <- g +
scale_fill_manual(values = colors, labels = names_colors, drop = FALSE, limits = names_colors, name = group.by,
guide = guide_legend(override.aes = list(alpha = 1), order = 2)) +
scale_color_manual(values = colors, labels = names_colors, drop = FALSE, limits = names_colors, name = group.by,
guide = guide_legend(override.aes = list(alpha = 1), order = 2))
} else {
stop("Only ROIs, spots, cells, molecules and tiles can be visualized with vrSpatialPlot!")
}
return(g)
}
####
## Spatial Feature Plot ####
####
#' vrSpatialFeaturePlot
#'
#' Plotting single/multiple features of spatially resolved cells, spots, and ROI on associated images from multiple assays in a VoltRon object.
#'
#' @param object a VoltRon object
#' @param features a set of features to be visualized, either from \link{vrFeatures} of raw or normalized data or columns of the \link{Metadata}.
#' @param combine.features whether to combine all features in one plot
#' @param group.by a column of metadata from \link{Metadata} used as grouping label for the spatial entities
#' @param plot.segments plot segments from \link{vrSegments} instead of points
#' @param n.tile should points be aggregated into tiles before visualization (see \link{geom_tile}). Applicable only for cells and molecules
#' @param norm if TRUE, the normalized data is used
#' @param log if TRUE, data features (excluding metadata features) will be log transformed
#' @param assay assay name (exp: Assay1) or assay class (exp: Visium, Xenium), see \link{SampleMetadata}.
#' if NULL, the default assay will be used, see \link{vrMainAssay}.
#' @param graph.name if not NULL, the spatial graph is with name \code{graph.name} is visualized as well, see \link{vrGraphNames}
#' @param ncol column wise number of plots, for \link{ggarrange}
#' @param nrow row wise number of plots, for \link{ggarrange}
#' @param font.size font size
#' @param pt.size point size
#' @param cell.shape the shape of the points representing cells, see \link{geom_point}
#' @param title.size title size of legend and plot
#' @param alpha alpha level of colors of visualized points and segments
#' @param keep.scale whether unify all scales for all features or not
#' @param label if TRUE, labels of ROIs will be visualized too
#' @param spatial the name of the main spatial system
#' @param channel the name of the channel associated with the image
#' @param background.color the color of plot background if a channel is not specified, or the spatial coord system doesnt have an image.
#' @param background (DEPRECATED) the background of the plot. Either an image name, see \link{vrImageNames} or a vector of length two with image name
#' and a channel name, see \link{vrImageChannelNames}. Type "black" or "white" for black or white backgrounds. if NULL, the main image (\link{vrMainSpatial})
#' and main channel (\link{vrMainChannel}) will be in the background. Otherwise the background will be grey.
#' @param reg TRUE if registered coordinates of the main image (\link{vrMainSpatial}) is requested
#' @param crop whether to crop an image of a spot assay to the extend of spots
#' @param scale.image if TRUE, background image will be scaled down to a low resolution (width: 1000px)
#' @param common.legend whether to use a common legend for all plots, see \link{ggarrange}
#' @param legend.loc the location of the legend, default is "right"
#' @param collapse.plots whether to combine all ggplots
#'
#' @import ggplot2
#' @importFrom ggpubr ggarrange
#'
#' @export
#'
vrSpatialFeaturePlot <- function(object, features, combine.features = FALSE, group.by = "label", plot.segments = FALSE, n.tile = 0,
norm = TRUE, log = FALSE, assay = NULL, graph.name = NULL, ncol = 2, nrow = NULL,
font.size = 2, pt.size = 2, cell.shape = 16, title.size = 10, alpha = 0.6, keep.scale = "feature", label = FALSE,
spatial = NULL, channel = NULL, background.color = NULL, background = NULL, reg = FALSE,
crop = FALSE, scale.image = TRUE, common.legend = FALSE, legend.loc = "right", collapse.plots = TRUE) {
# check object
if(!inherits(object, "VoltRon"))
stop("Please provide a VoltRon object!")
# sample metadata
sample.metadata <- SampleMetadata(object)
# get assay names
assay_names <- vrAssayNames(object, assay = assay)
# get entity type and metadata
metadata <- Metadata(object, assay = assay)
# calculate limits for plotting, all for making one scale, feature for making multiple
limits <- Map(function(feat){
range_feat <- Map(function(assy){
spatialpoints <- vrSpatialPoints(object[[assy]])
if(feat %in% vrFeatures(object, assay = assy)){
data <- vrData(object[[assy]], features = feat, norm = norm)
if(log)
data <- log1p(data)
# return(getRange(data, na.rm = TRUE, finite = TRUE))
return(getRange(data, na.rm = TRUE))
} else {
metadata <- Metadata(object, assay = assy)
if(feat %in% colnames(metadata)){
# return(getRange(metadata[,feat], na.rm = TRUE, finite = TRUE))
return(getRange(metadata[,feat], na.rm = TRUE))
} else {
stop("Feature '", feat, "' cannot be found in data or metadata!")
}
}
}, assay_names)
if(keep.scale == "all"){
range_feat_all <- c(do.call(min, range_feat), do.call(max, range_feat))
range_feat <- Map(function(assy) return(range_feat_all), assay_names)
}
if(!keep.scale %in% c("all", "feature")){
stop("keep.scale should be either 'all' or 'feature', check help(vrSpatialFeaturePlot)")
}
return(range_feat)
}, features)
# configure titles
assay_title <- as.list(apply(sample.metadata[assay_names,], 1, function(x) paste(x["Sample"], x["Layer"], sep = ", ")))
names(assay_title) <- assay_names
feature_title <- as.list(features)
names(feature_title) <- features
plot_title <- assay_title
legend_title <- feature_title
# for each assay
i <- 1
gg <- list()
for(assy in assay_names){
# get assay
cur_assay <- object[[assy]]
cur_metadata <- subset_metadata(metadata, assays = assy)
p_title <- plot_title[[assy]]
# get graph
if(!is.null(graph.name)){
if(graph.name %in% vrGraphNames(object)){
graph <- vrGraph(object, assay = assy, graph.type = graph.name)
} else {
stop("the graph with named '", graph.name, "' was not found in the graph list!")
}
} else {
graph <- NULL
}
# make one plot by combining features
if(combine.features){
if(vrAssayTypes(cur_assay) != "cell"){
stop("Currently, combining features can only be used for cell based assays")
} else {
l_title <- ""
gg[[i]] <- vrSpatialFeaturePlotCombined(assay = cur_assay, metadata = cur_metadata, features = features, plot.segments = plot.segments,
n.tile = n.tile, graph = graph, group.by = group.by, norm = norm, log = log,
font.size = font.size, pt.size = pt.size, title.size = title.size, alpha = alpha, cell.shape = cell.shape,
label = label, plot_title = p_title, legend_title = l_title, background = background, reg = reg, crop = crop)
}
# make individual feature plots
} else {
# for each feature
for(feat in features){
l_title <- legend_title[[feat]]
gg[[i]] <- vrSpatialFeaturePlotSingle(assay = cur_assay, metadata = cur_metadata, feature = feat, plot.segments = plot.segments,
n.tile = n.tile, graph = graph, limits = limits[[feat]][[assy]], group.by = group.by, norm = norm, log = log,
font.size = font.size, pt.size = pt.size, title.size = title.size, alpha = alpha, cell.shape = cell.shape,
label = label, plot_title = p_title, legend_title = l_title, spatial = spatial, channel = channel,
background.color = background.color, background = background, reg = reg, crop = crop)
i <- i + 1
}
}
}
if(collapse.plots){
# return a list of plots or a single one
if(length(features) > 1 && length(assay_names) > 1){
return(suppressWarnings(ggpubr::ggarrange(plotlist = gg, ncol = length(features), nrow = length(assay_names))))
} else if(length(features) > 1 && length(assay_names) == 1){
if(length(gg) < ncol) ncol <- length(gg)
return(suppressWarnings(ggpubr::ggarrange(plotlist = gg, ncol = ncol, nrow = ceiling(length(gg)/ncol))))
} else if(length(features) == 1 && length(assay_names) > 1){
if(length(gg) < ncol) ncol <- length(gg)
return(suppressWarnings(ggpubr::ggarrange(plotlist = gg, ncol = ncol, nrow = ceiling(length(gg)/ncol),
common.legend = common.legend, legend = legend.loc)))
} else {
return(suppressWarnings(gg[[1]] + theme(legend.position=legend.loc)))
}
} else {
for(i in seq_len(length(features)*length(assay_names))){
gg[[i]] <- gg[[i]] + theme(legend.position=legend.loc)
}
return(suppressWarnings(gg))
}
}
#' vrSpatialFeaturePlotSingle
#'
#' A single Spatial Feature plot of VoltRon object
#'
#' @param assay vrAssay object
#' @param metadata the metadata associated with the assay
#' @param feature a feature, either from the rows of rawdata, normdata or columns of the metadata
#' @param plot.segments plot segments from \link{vrSegments} instead of points
#' @param n.tile should points be aggregated into tiles before visualization (see \link{geom_tile}). Applicable only for cells and molecules
#' @param graph if not NULL, the graph is added to the plot
#' @param limits limits of the legend of the plot
#' @param group.by a column of metadata from \link{Metadata} used as grouping label for the spatial entities
#' @param norm if TRUE, the normalized data is used
#' @param log if TRUE, data features (excluding metadata features) will be log transformed
#' @param font.size font sizes
#' @param pt.size point size
#' @param cell.shape the shape of the points representing cells, see \link{geom_point}
#' @param title.size title size of legend and plot
#' @param alpha alpha level of colors of visualized points and segments
#' @param label if TRUE, labels of ROIs will be visualized too
#' @param plot_title the main title of the single plot
#' @param legend_title the legend title of the single plot
#' @param spatial the name of the main spatial system
#' @param channel the name of the channel associated with the image
#' @param background.color the color of plot background if a channel is not specified, or the spatial coord system doesnt have an image.
#' @param background (DEPRECATED) the background of the plot. Either an image name, see \link{vrImageNames} or a vector of length two with image name
#' and a channel name, see \link{vrImageChannelNames}. Type "black" or "white" for black or white backgrounds. if NULL, the main image (\link{vrMainSpatial})
#' and main channel (\link{vrMainChannel}) will be in the background. Otherwise the background will be grey.
#' @param reg TRUE if registered coordinates of the main image (\link{vrMainSpatial}) is requested
#' @param crop whether to crop an image of a spot assay to the extend of spots
#' @param scale.image if TRUE, background image will be scaled down to a low resolution (width: 1000px)
#'
#' @import ggplot2
#' @importFrom ggrepel geom_label_repel
#' @importFrom igraph get.data.frame
#' @importFrom dplyr arrange
#'
#' @noRd
vrSpatialFeaturePlotSingle <- function(assay, metadata, feature, plot.segments = FALSE, n.tile = 0, graph = NULL, limits,
group.by = "label", norm = TRUE, log = FALSE, font.size = 2, pt.size = 2, cell.shape = 16,
title.size = 10, alpha = 0.6, label = FALSE, plot_title = NULL, legend_title = NULL,
spatial = NULL, channel = NULL, background.color = NULL, background = NULL, reg = FALSE,
crop = FALSE, scale.image = TRUE){
# get image information and plotting features
midpoint <- sum(limits)/2
# get data
data_features <- feature[feature %in% vrFeatures(assay)]
if(length(data_features) > 0){
normdata <- vrData(assay, features = feature, norm = norm)
if(log)
normdata <- log1p(normdata)
}
if(feature %in% data_features){
if(inherits(normdata, "IterableMatrix")){
cur_score <- as.matrix(normdata[feature,])[1,]
} else {
cur_score <- normdata[feature,]
names(cur_score) <- colnames(normdata)
}
} else {
cur_score <- metadata[,feature]
if(!is.null(rownames(metadata))){
names(cur_score) <- rownames(metadata)
} else {
names(cur_score) <- as.vector(metadata$id)
}
}
# plot
g <- ggplot()
# add image and background
image <- vrSpatialPlotImage(g, assay, background, scale.image, spatial = spatial,
channel = channel, background.color = background.color)
g <- image$plot
info <- image$info
background.color <- image$background.color
spatial_name <- image$spatial
scale_factors <- image$scale_factors
# coords
coords <- vrCoordinates(assay, spatial_name = spatial_name, reg = reg)
if(!inherits(coords, "IterableMatrix")){
coords <- as.data.frame(coords)
}
coords <- coords/scale_factors
segments <- vrSegments(assay, spatial_name = spatial_name)
# add score
coords <- as.matrix(coords)
coords <- data.frame(coords, cur_score[rownames(coords), drop = FALSE])
colnames(coords)[length(colnames(coords))] <- "score"
# set up the limits
g <- vrSpatialExtent(g, assay, coords, crop, info)
# add points or segments
if(vrAssayTypes(assay) == "ROI" && !is.null(segments)){
polygon_data <- NULL
circle_data <- NULL
for(i in seq_len(length(segments))){
cur_data <- as.data.frame(cbind(segments[[i]], names(segments)[i], coords$score[i]))
if(nrow(segments[[i]]) > 1){
colnames(cur_data) <- c("id", "x", "y", "segment", "score")
cur_data[,c("x", "y")] <- cur_data[,c("x", "y")]/scale_factors
polygon_data <- as.data.frame(rbind(polygon_data, cur_data))
} else {
colnames(cur_data) <- c("id", "x", "y", "rx", "ry", "segment", "score")
cur_data[,c("x", "y","rx", "ry")] <- cur_data[,c("x", "y","rx", "ry")]/scale_factors
circle_data <- as.data.frame(rbind(circle_data, cur_data))
}
}
if(!is.null(geom_polygon)){
g <- g +
geom_polygon(aes(x = x, y = y, fill = score, group = segment), data = polygon_data, alpha = alpha)
}
if(!is.null(circle_data)){
if(!requireNamespace('ggforce'))
stop("Please install ggforce package!: install.packages('ggforce')")
g <- g +
ggforce::geom_ellipse(aes(x0 = as.numeric(x), y0 = as.numeric(y), a = as.numeric(rx), b = as.numeric(ry), angle = 0,
fill = score, group = segment), data = circle_data, lwd = 0, alpha = alpha)
}
g <- g +
scale_fill_gradientn(name = legend_title,
colors=c("dodgerblue2", "white", "yellow3"),
values=rescale_numeric(c(limits[1], midpoint, limits[2])), limits = limits)
} else if(vrAssayTypes(assay) == "spot"){
# rasterize if requested or needed
if(n.tile > 0 || nrow(coords) > 50000){
if(n.tile == 0)
n.tile <- 1000
g <- vrFeaturePlotTiling(g = g, data = coords, legend_title = legend_title, n.tile = n.tile, alpha = alpha, type = "spatial")
} else {
spot.type <- vrAssayParams(assay, param = "spot.type")
spot.type <- ifelse(is.null(spot.type), "circle", spot.type)
g <- g +
geom_spot(mapping = aes(x = x, y = y, fill = score), coords, shape = 21, alpha = alpha,
spot.radius = vrAssayParams(assay, param = "vis.spot.radius")/scale_factors,
spot.type = spot.type) +
scale_fill_gradientn(name = legend_title,
colors=c("dodgerblue3", "yellow", "red"),
values=rescale_numeric(c(limits[1], midpoint, limits[2])), limits = limits)
}
} else if(vrAssayTypes(assay) %in% c("cell", "tile")) {
if(plot.segments){
if(length(segments) == 0) {
stop("No Segments are available in this assay!")
} else {
polygon_data <- do.call(rbind,segments)
polygon_data[,c("x", "y")] <- polygon_data[,c("x", "y")]/scale_factors
# len_segments <- sapply(segments, nrow, simplify = TRUE)
len_segments <- vapply(segments, nrow, numeric(1))
polygon_data <- data.frame(polygon_data, segment = rep(names(segments), len_segments), score = rep(coords$score, len_segments))
g <- g +
geom_polygon(aes(x = x, y = y, fill = score, group = segment), data = polygon_data, alpha = alpha)
}
g <- g +
scale_fill_gradientn(name = legend_title,
colors=c("dodgerblue2", "white", "yellow3"),
values=rescale_numeric(c(limits[1], midpoint, limits[2])), limits = limits)
} else {
# rasterize if requested or needed
if(n.tile > 0 || nrow(coords) > 50000){
if(n.tile == 0)
n.tile <- 1000
g <- vrFeaturePlotTiling(g = g, data = coords, legend_title = legend_title, n.tile = n.tile, alpha = alpha, type = "spatial")
} else {
g <- g +
geom_point(mapping = aes(x = x, y = y, fill = score, color = score), dplyr::arrange(coords, score), shape = cell.shape, size = rel(pt.size), alpha = alpha) +
scale_colour_gradientn(name = legend_title,
colors=c("dodgerblue2", "white", "yellow3"),
values=rescale_numeric(c(limits[1], midpoint, limits[2])), limits = limits, aesthetics = c("fill", "colour"))
}
# add if a graph exists
if(!is.null(graph)){
graph.df <- igraph::get.data.frame(graph)
graph.df$from.x <- coords$x[match(graph.df$from, rownames(coords))]
graph.df$from.y <- coords$y[match(graph.df$from, rownames(coords))]
graph.df$to.x <- coords$x[match(graph.df$to, rownames(coords))]
graph.df$to.y <- coords$y[match(graph.df$to, rownames(coords))]
g <- g +
geom_segment(data = graph.df, mapping = aes(x=from.x,xend = to.x, y=from.y,yend = to.y), alpha = 0.5, color = ifelse(background == "black", "grey", "black"))
}
}
} else {
stop("Feature plotting is not possible for molecule assays!")
}
# more visualization parameters
g <- g +
ggtitle(plot_title) + theme(plot.title = element_text(hjust = 0.5, margin=margin(0,0,0,0), size = title.size),
panel.grid.minor = element_blank(), panel.grid.major = element_blank(),
axis.line=element_blank(),axis.text.x=element_blank(), axis.text.y=element_blank(),
axis.ticks=element_blank(), axis.title.x=element_blank(), axis.title.y=element_blank(),
legend.key.size = unit(title.size, "points"), legend.title = element_text(size=title.size),
legend.margin = margin(0,0,0,0))
# visualize labels
if(label && vrAssayTypes(assay) == "ROI"){
if(group.by %in% colnames(metadata)){
coords[[group.by]] <- as.vector(metadata[,group.by])
} else {
stop("The column ", group.by, " was not found in the metadata!")
}
g <- g + ggrepel::geom_label_repel(mapping = aes(x = .data[["x"]], y = .data[["y"]], label = .data[[group.by]]), coords,
box.padding = 0.5, size = font.size, direction = "both", seed = 1)
}
# return data
return(g)
}
#' vrSpatialFeaturePlotSingle
#'
#' A single Spatial Feature plot of VoltRon object
#'
#' @param assay vrAssay object
#' @param metadata the metadata associated with the assay
#' @param features features, either from the rows of rawdata, normdata or columns of the metadata
#' @param plot.segments plot segments from \link{vrSegments} instead of points
#' @param n.tile should points be aggregated into tiles before visualization (see \link{geom_tile}). Applicable only for cells and molecules
#' @param graph if not NULL, the graph is added to the plot
#' @param limits limits of the legend of the plot
#' @param group.by a column of metadata from \link{Metadata} used as grouping label for the spatial entities
#' @param norm if TRUE, the normalized data is used
#' @param log if TRUE, data features (excluding metadata features) will be log transformed
#' @param font.size font sizes
#' @param pt.size point size
#' @param cell.shape the shape of the points representing cells, see \link{geom_point}
#' @param title.size title size of legend and plot
#' @param alpha alpha level of colors of visualized points and segments
#' @param label if TRUE, labels of ROIs will be visualized too
#' @param plot_title the main title of the single plot
#' @param legend_title the legend title of the single plot
#' and a channel name, see \link{vrImageChannelNames}. Type "black" or "white" for black or white backgrounds. if NULL, the main image (\link{vrMainSpatial})
#' and main channel (\link{vrMainChannel}) will be in the background. Otherwise the background will be grey.
#' @param background the background of the plot. Either an image name, see \link{vrImageNames} or a vector of length two with image name
#' and a channel name, see \link{vrImageChannelNames}. Type "black" or "white" for black or white backgrounds. if NULL, the main image (\link{vrMainSpatial})
#' and main channel (\link{vrMainChannel}) will be in the background. Otherwise the background will be grey.
#' @param reg TRUE if registered coordinates of the main image (\link{vrMainSpatial}) is requested
#' @param crop whether to crop an image of a spot assay to the extend of spots
#' @param scale.image if TRUE, background image will be scaled down to a low resolution (width: 1000px)
#'
#' @import ggplot2
#' @importFrom ggrepel geom_label_repel
#' @importFrom igraph get.data.frame
#' @importFrom dplyr arrange
#'
#' @noRd
vrSpatialFeaturePlotCombined <- function(assay, metadata, features, plot.segments = FALSE, n.tile = 0, graph = NULL, group.by = "label", norm = TRUE, log = FALSE,
font.size = 2, pt.size = 2, cell.shape = 16, title.size = 10, alpha = 0.6, label = FALSE, plot_title = NULL,
legend_title = NULL, background = NULL, reg = FALSE, crop = FALSE, scale.image = TRUE){
# plot
g <- ggplot()
scale_factors <- 1
# add image
info <- NULL
spatial_name <- vrMainSpatial(assay)
# data
coords <- as.data.frame(vrCoordinates(assay, spatial_name = spatial_name, reg = reg))
coords <- coords/scale_factors
segments <- vrSegments(assay, spatial_name = spatial_name)
data_features <- features[features %in% vrFeatures(assay)]
if(length(data_features) > 0){
normdata <- vrData(assay, features = data_features, norm = norm)
if(log)
normdata <- log1p(normdata)
}
# calculate limits for plotting, all for making one scale, feature for making multiple
limits <- Map(function(feat){
if(feat %in% vrFeatures(assay)){
return(range(normdata[feat, ]))
} else {
if(feat %in% colnames(metadata)){
return(range(metadata[,feat]))
} else {
stop("Feature '", feat, "' cannot be found in data or metadata!")
}
}
}, features)
# set up the limits
g <- vrSpatialExtent(g, assay, coords, crop, info)
# configure titles
feature_title <- as.list(features)
names(feature_title) <- features
legend_title <- feature_title
# for each feature
i <- 1
gg <- list()
all_data <- NULL
colors <- get_rasterization_colors(length(features))
for(feat in features){
# get data
if(feat %in% data_features){
coords$score <- normdata[feat,]
} else {
coords$score <- metadata[,feat]
}
# get image information and plotting features
midpoint <- sum(limits[[feat]])/2
# add points, rasterize if requested or needed
if(n.tile > 0 || nrow(coords) > 50000){
if(n.tile == 0)
n.tile <- 1000
g_single <- ggplot() +
stat_summary_2d(mapping = aes(x = x, y = y, z = score), fun = mean, data = coords, geom = "tile", bins = n.tile, drop = TRUE) +
scale_fill_gradientn(name = legend_title[[feat]],
colors=c("grey97", colors[i]),
values=rescale_numeric(c(limits[[feat]][1], limits[[feat]][2])), limits = limits[[feat]])
all_data <- rbind(all_data,
data.frame(layer_data(g_single), color_group = colors[i]))
} else {
g_single <- ggplot() +
geom_point(mapping = aes(x = x, y = y, color = score), coords, shape = 16, size = pt.size) +
scale_color_gradientn(name = legend_title[[feat]],
colors=c("grey97", colors[i]),
values=rescale_numeric(c(limits[[feat]][1], limits[[feat]][2])), limits = limits[[feat]])
all_data <- rbind(all_data,
data.frame(layer_data(g_single), value = coords$score, color_group = colors[i]))
}
# add graph to list
gg[[i]] <- g_single
i <- i + 1
}
# combine feature plots
g <- vrSpatialFeatureCombinePlot(g, all_data, n.tile, coords, features)
# add if a graph exists
if(!is.null(graph)){
graph.df <- igraph::get.data.frame(graph)
graph.df$from.x <- coords$x[match(graph.df$from, rownames(coords))]
graph.df$from.y <- coords$y[match(graph.df$from, rownames(coords))]
graph.df$to.x <- coords$x[match(graph.df$to, rownames(coords))]
graph.df$to.y <- coords$y[match(graph.df$to, rownames(coords))]
g <- g +
geom_segment(data = graph.df, mapping = aes(x=from.x,xend = to.x, y=from.y,yend = to.y), alpha = 0.5, color = ifelse(background == "black", "grey", "black"))
}
# more visualization parameters
g <- g +
ggtitle(plot_title) + theme(plot.title = element_text(hjust = 0.5, margin=margin(0,0,0,0), size = title.size),
panel.grid.minor = element_blank(), panel.grid.major = element_blank(),
axis.line=element_blank(),axis.text.x=element_blank(), axis.text.y=element_blank(),
axis.ticks=element_blank(), axis.title.x=element_blank(), axis.title.y=element_blank(),
legend.key.size = unit(title.size, "points"), legend.title = element_text(size=title.size),
legend.margin = margin(0,0,0,0))
# background
g <- g +
theme(panel.background = element_rect(fill = "grey97", colour = "grey97", linewidth = 0.5, linetype = "solid"))
# return plot
return(g)
}
#' vrSpatialFeatureCombinePlot
#'
#' A different plotting setting for vrSpatialFeaturePlot when combine.feature=TRUE
#'
#' @param g ggplot object
#' @param all_data summary data
#' @param n.tile should points be aggregated into tiles before visualization (see \link{geom_tile}). Applicable only for cells and molecules
#' @param datax original plotting data
#' @param features features
#'
#' @import ggplot2
#'
#' @noRd
vrSpatialFeatureCombinePlot <- function(g, all_data, n.tile, coords, features){
# tiling or not
if(n.tile > 0 || nrow(coords) > 50000){
if(n.tile == 0)
n.tile <- 1000
all_data <- all_data %>% group_by(x,y) %>%
summarize(fill = fill[which.max(value)],
group = color_group[which.max(value)],
value = value[which.max(value)])
key_table <- all_data[,c("fill", "group", "value")] %>%
dplyr::group_by(group) %>%
dplyr::summarise(fill = fill[which.max(value)], value = max(value))
g.combined <- g +
ggplot2::geom_tile(data = as.data.frame(all_data), aes(x = x, y = y, fill = fill)) +
ggplot2::scale_fill_identity("", labels = features, breaks = key_table$fill, guide = "legend")
} else {
all_data <- all_data %>% group_by(x,y) %>%
dplyr::summarise(color = colour[which.max(value)],
group = color_group[which.max(value)],
value = value[which.max(value)])
key_table <- all_data[,c("color", "group", "value")] %>%
dplyr::group_by(group) %>%
dplyr::summarise(color = color[which.max(value)], value = max(value))
g.combined <- g +
ggplot2::geom_point(data = as.data.frame(all_data), aes(x = x, y = y, color = color)) +
ggplot2::scale_color_identity("", labels = features, breaks = key_table$color, guide = "legend")
}
g.combined
}
####
## Spatial Auxiliary ####
####
#' vrSpatialPlotImage
#'
#' setting background with or without the image
#'
#' @param g a ggplot object
#' @param assay vrAssay object
#' @param spatial the name of the main spatial system
#' @param channel the name of the channel associated with the image
#' @param background.color the color of plot background if a channel is not specified, or the spatial coord system doesnt have an image.
#' @param background (DEPRECATED) the background of the plot. Either an image name, see \link{vrImageNames} or a vector of length two with image name
#' and a channel name, see \link{vrImageChannelNames}. Type "black" or "white" for black or white backgrounds. if NULL, the main image (\link{vrMainSpatial})
#' and main channel (\link{vrMainChannel}) will be in the background. Otherwise the background will be grey.
#' @param scale.image if TRUE, background image will be scaled down to a low resolution (width: 1000px)
#'
#' @import ggplot2
#' @importFrom magick image_read
#'
#' @noRd
vrSpatialPlotImage <- function(g, assay, background, scale.image, spatial = NULL, channel = NULL, background.color = NULL){
# parameters
scale_factors <- 1
# check background, spatial and channel
if(!is.null(background)){
warning("The use of 'background' is deprecated, please use 'spatial' to call spatial coordinate systems and 'channel' to get channels of associated images.")
warning("Please use 'background.color' to draw spatial plots with custom colors in the background.")
if(length(background) == 2) {
channel <- background[2]
} else {
channel <- NULL
}
spatial <- background[1]
} else {
if(is.null(spatial)){
spatial <- vrMainSpatial(assay)
}
}
# get image and adjust size
if(spatial %in% vrSpatialNames(assay) && is.null(background.color)){
# get image
image <- suppressWarnings({vrImages(assay, name = spatial, channel = channel, as.raster = TRUE)})
if(!is.null(image) & !inherits(image, "Image_Array")){
image <- magick::image_read(image)
}
if(!is.null(image)){
info <- getImageInfo(image)
if(info$width > 1000 && scale.image){
image <- resize_Image(image, geometry = "1000x")
scale_factors <- info$width/1000
info <- getImageInfo(image)
}
g <- g +
ggplot2::annotation_raster(image, 0, info$width, info$height, 0, interpolate = FALSE)
} else {
info <- NULL
}
} else {
info <- NULL
}
# background color
if(!is.null(background.color)){
g <- g +
theme(panel.background = element_rect(fill = background.color, colour = background.color, linewidth = 0.5, linetype = "solid"))
} else{
if(is.null(info)){
g <- g +
theme(panel.background = element_rect(fill = "grey97", colour = "grey97", linewidth = 0.5, linetype = "solid"))
} else{
g <- g +
theme(panel.background = element_blank())
}
}
list(plot = g, info = info, background.color = background.color, spatial = spatial, scale_factors = scale_factors)
}
#' vrSpatialExtent
#'
#' setting extent for spatial plots
#'
#' @param g a ggplot object
#' @param assay vrAssay object
#' @param coords coordinates
#' @param crop whether to crop an image of a spot assay to the extend of spots
#' @param image info
#'
#' @importFrom ggplot2 xlim ylim coord_fixed
#'
#' @noRd
vrSpatialExtent <- function(g, assay, coords, crop, info){
if(vrAssayTypes(assay) == "spot"){
if(crop){
g <- g +
ggplot2::coord_fixed(xlim = range(coords$x), ylim = range(coords$y))
} else {
if(!is.null(info)){
g <- g +
ggplot2::coord_fixed(xlim = c(0,info$width), ylim = c(0,info$height))
}
}
} else {
if(crop){
g <- g +
ggplot2::coord_fixed(xlim = range(coords$x), ylim = range(coords$y))
} else {
if(!is.null(info)){
g <- g +
ggplot2::xlim(0,info$width) + ggplot2::ylim(0, info$height)
}
}
}
g
}
#' @import ggplot2
#' @importFrom grid pointsGrob unit gpar
#' @importFrom rlang list2
#'
#' @noRd
NULL
#' geom_spot
#'
#' @import ggplot2
#' @importFrom grid pointsGrob unit gpar
#' @importFrom rlang list2
#'
#' @noRd
geom_spot <- function (mapping = NULL, data = NULL, stat = "identity", position = "identity",
..., na.rm = FALSE, show.legend = NA, inherit.aes = TRUE) {
ggplot2::layer(data = data, mapping = mapping, stat = stat, geom = GeomSpot,
position = position, show.legend = show.legend, inherit.aes = inherit.aes,
params = rlang::list2(na.rm = na.rm, ...))
}
#' GeomSpot (ggproto class)
#'
#' @import ggplot2
#' @importFrom grid pointsGrob unit gpar
#'
#' @noRd
GeomSpot <- ggplot2::ggproto("GeomSpot",
Geom,
required_aes = c("x", "y"),
non_missing_aes = c("size", "shape", "colour"),
default_aes = aes(
shape = 21,
colour = "black",
size = 1.5,
fill = NA,
spot.radius = 1,
spot.type = "circle",
alpha = NA,
stroke = 0.5
),
draw_panel = function(self, data, panel_params, coord, na.rm = FALSE) {
if (is.character(data$shape)) {
data$shape <- translate_shape_string(data$shape)
}
if(!is.null(data$spot.type)){
if(all(data$spot.type == "square")){
data$shape <- 15
data$spot.radius <- data$spot.radius/2
}
}
coords <- coord$transform(data, panel_params)
stroke_size <- coords$stroke
stroke_size[is.na(stroke_size)] <- 0
xrange <- panel_params$x.range[2] - panel_params$x.range[1]
yrange <- panel_params$y.range[2] - panel_params$y.range[1]
mainrange <- min(xrange, yrange)
spot.radius <- data$spot.radius/mainrange
ggname("geom_spot",
grid::pointsGrob(
coords$x, coords$y,
pch = coords$shape,
size = grid::unit(spot.radius, "npc"),
gp = grid::gpar(
col = ggplot2::alpha(coords$colour, coords$alpha),
fill = ggplot2::alpha(coords$fill, coords$alpha),
fontsize = coords$size * .pt + stroke_size * .stroke / 2,
lwd = coords$stroke * .stroke / 2
)
)
)
},
draw_key = draw_key_point
)
#' addGraph
#'
#' add graph to a spatial plot
#'
#' @param coords coordinates data
#' @param graph if not NULL, the graph is added to the plot
#' @param background the background of the plot. Either an image name, see \link{vrImageNames} or a vector of length two with image name
#'
#' @import ggplot2
#' @importFrom igraph as_data_frame
#'
#' @noRd
addGraph <- function(graph, coords, background){
# graph.df <- igraph::get.data.frame(graph)
graph.df <- igraph::as_data_frame(graph)
graph.df$from.x <- coords$x[match(graph.df$from, rownames(coords))]
graph.df$from.y <- coords$y[match(graph.df$from, rownames(coords))]
graph.df$to.x <- coords$x[match(graph.df$to, rownames(coords))]
graph.df$to.y <- coords$y[match(graph.df$to, rownames(coords))]
geom_segment(data = graph.df, mapping = aes(x=from.x,xend = to.x, y=from.y,yend = to.y),
alpha = 0.5, color = background)
}
####
## Spatial Neighborhood Plot ####
####
#' vrNeighbourhoodEnrichmentPlot
#'
#' Plotting results of \link{vrNeighbourhoodEnrichment}
#'
#' @param results The results from \link{vrNeighbourhoodEnrichment}
#' @param assay assay name (exp: Assay1), see \link{SampleMetadata}.
#' @param type the type of spatial test. Either "assoc" for association test or "segreg" for segregation test
#'
#' @export
#'
vrNeighbourhoodEnrichmentPlot <- function(results, assay = NULL, type = c("assoc", "segreg")){
# check Seurat package
if(!requireNamespace('circlize'))
stop("Please install circlize package for coloring the Heatmap!: install.packages('circlize')")
# get assay names
if(!is.null(assay)){
if(length(assay) > 1){
stop("Spatial enrichment plots of only one Assay can be visualized at a time")
}
} else {
stop("assay should be defined, e.g. assay == 'Assay1'")
}
results <- results[results$AssayID == assay,]
# get type
if(length(type) > 1)
type <- "assoc"
# make vis matrices
results$p_adj <- ifelse(results$value > 0, results$p_assoc, results$p_segreg)
results$cell1 <- results$from_value
results$cell2 <- results$to_value
entities <- unique(c(results$cell1,results$cell2))
mat <- matrix(nrow = length(entities), ncol = length(entities))
rownames(mat) <- colnames(mat) <- entities
mat_padj <- mat
for(i in seq_len(nrow(results))){
mat_padj[which(results$cell1[i]==entities), which(results$cell2[i]==entities)] <-
mat_padj[which(results$cell2[i]==entities), which(results$cell1[i]==entities)] <- ifelse(type == "assoc", results$p_assoc_adj[i], results$p_segreg_adj[i])
mat[which(results$cell1[i]==entities), which(results$cell2[i]==entities)] <-
mat[which(results$cell2[i]==entities), which(results$cell1[i]==entities)] <- ifelse(type == "assoc", results$value[i], -1*results$value[i])
}
# color fun
limits <- c(0,round(quantile(mat[mat > 0], probs = 0.99),2))
enrich_col_fun = circlize::colorRamp2(limits, c("white", "red"))
# make heatmap
ht <- ComplexHeatmap::Heatmap(mat_padj,
cluster_rows = TRUE, show_row_dend = FALSE,
cluster_columns = TRUE, show_column_dend = TRUE,
column_names_rot = 45,
col = enrich_col_fun,
heatmap_legend_param = list(
at = limits,
title = "Score"),
cell_fun = function(j, i, x, y, width, height, fill) {
grid::grid.rect(x = x, y = y, width = width, height = height,
gp = gpar(col = "black", fill = "gray75"))
grid::grid.circle(x = x, y = y, r = (1-mat_padj[i, j])/2 * (min(grid::unit.c(width, height))*0.90),
gp = gpar(fill = enrich_col_fun(mat[i, j]), col = "black"))
})
# extra heatmap legend
lgd <- ComplexHeatmap::Legend(labels = c("1.00", "0.75", "0.50", "0.25"), title = "p.adj", type = "points", pch = 16,
size = unit(seq_len(5), 'mm'), direction = 'vertical', legend_gp = gpar(col = "black"), background = 'white')
# draw heatmap
ComplexHeatmap::draw(ht, annotation_legend_list = lgd)
}
####
# Embeddings plots ####
####
####
## Embedding Identity Plot ####
####
#' vrEmbeddingPlot
#'
#' Plotting embeddings of cells and spots on associated images from multiple assays in a VoltRon object.
#'
#' @param object a VoltRon object
#' @param embedding the embedding type, i.e. pca, umap etc.
#' @param group.by a column of metadata from \link{Metadata} used as grouping label for the spatial entities
#' @param group.ids a subset of categories defined in metadata column from \code{group.by}
#' @param split.by a column of metadata from \link{Metadata} used as splitting spatial entities into ggplot panels, see \link{facet_wrap}
#' @param colors the color set for group.by. Should be of the same size of group.id (if specified) or unique elements in group.by
#' @param n.tile should points be aggregated into tiles before visualization (see \link{geom_tile}). Applicable only for cells and molecules
#' @param assay assay name (exp: Assay1) or assay class (exp: Visium, Xenium), see \link{SampleMetadata}.
#' if NULL, the default assay will be used, see \link{vrMainAssay}.
#' @param ncol column wise number of plots, for \link{facet_wrap}
#' @param nrow row wise number of plots, for \link{facet_wrap}
#' @param font.size font size of labels, if label is TRUE
#' @param pt.size point size
#' @param label if TRUE, the labels of the ROI assays will be visualized
#' @param common.legend whether to use a common legend for all plots, see \link{ggarrange}
#' @param collapse.plots whether to combine all ggplots
#'
#' @import ggplot2
#' @importFrom stats aggregate
#'
#' @export
#'
vrEmbeddingPlot <- function(object, embedding = "pca", group.by = "Sample", group.ids = NULL, split.by = NULL, colors = NULL, n.tile = 0, assay = NULL,
ncol = 2, nrow = NULL, font.size = 5, pt.size = 1, label = FALSE, common.legend = TRUE, collapse.plots = TRUE) {
# check object
if(!inherits(object, "VoltRon"))
stop("Please provide a VoltRon object!")
# sample metadata
sample.metadata <- SampleMetadata(object)
# get assay names
assay_names <- vrAssayNames(object, assay = assay)
# get entity type and metadata
metadata <- Metadata(object, assay = assay)
# grep assays from metadata
metadata <- subset_metadata(metadata, assays = assay_names)
# check group.by
if(!group.by %in% colnames(metadata)){
stop("Column ", group.by, " is not found in metadata!")
}
# adjust group.ids
if(is.null(group.ids)){
group.ids <- unique(metadata[[group.by]])
}
# check group.id
levels_group.ids <- as.character(group.ids)
if(all(!is.na(suppressWarnings(as.numeric(levels_group.ids))))){
levels_group.ids <- sort(as.numeric(levels_group.ids))
}
group.ids <- factor(group.ids, levels = levels_group.ids)
# check colors
if(!is.null(colors)){
if(length(colors) != length(levels_group.ids)){
stop("Please provide colors whose length is equal to group.ids (if specified) or unique elements of group.by")
}
if(is.null(names(colors))){
stop("Name of each color has to be specified given as names(colors)")
}
} else{
colors <- hue_pal(length(group.ids))
names(colors) <- group.ids
}
# plotting features
datax <- data.frame(vrEmbeddings(object, assay = assay_names, type = embedding))
datax <- datax[,seq_len(2)]
colnames(datax) <- c("x", "y")
if(group.by %in% colnames(metadata)){
if(inherits(metadata, "data.table")){
datax[[group.by]] <- metadata[,get(names(metadata)[which(colnames(metadata) == group.by)])]
} else {
if(!is.null(rownames(metadata))){
datax[[group.by]] <- as.factor(metadata[rownames(datax),group.by])
} else{
datax[[group.by]] <- as.factor(as.vector(metadata[match(rownames(datax), as.vector(metadata$id)),group.by]))
}
}
} else {
stop("Column ", group.by, " cannot be found in metadata!")
}
# subset group.by using group.id
datax <- droplevels(datax[datax[[group.by]] %in% group.ids,])
datax[[group.by]] <- factor(datax[[group.by]], levels = group.ids)
# get split.by
if(!is.null(split.by)){
if(split.by %in% colnames(metadata)){
if(inherits(metadata, "data.table")){
datax[[split.by]] <- metadata[,get(names(metadata)[which(colnames(metadata) == split.by)])]
} else {
datax[[split.by]] <- as.factor(metadata[rownames(datax),split.by])
}
} else {
stop("Column ", split.by, " cannot be found in metadata!")
}
}
# plot
g <- ggplot()
# add points, rasterize if requested or needed
if(n.tile > 0 || nrow(datax) > 50000){
if(n.tile == 0)
n.tile <- 1000
g <- vrGroupPlotTiling(g = g, data = datax, group.by = group.by, n.tile = n.tile)
} else {
# g <- g +
# geom_point(mapping = aes_string(x = "x", y = "y", color = group.by), datax, shape = 16, size = pt.size)
g <- g +
geom_point(mapping = aes(x = .data[["x"]], y = .data[["y"]],
color = .data[[group.by]], fill = .data[[group.by]]), datax, shape = 16, size = pt.size)
}
g <- g +
scale_color_manual(values = colors, labels = names(colors), drop = FALSE, limits = names(colors)) +
scale_fill_manual(values = colors, labels = names(colors), drop = FALSE, limits = names(colors)) +
guides(color = guide_legend(override.aes=list(size = 2)))
# more visualization parameters
g <- g +
theme_classic() +
xlab(paste0(toupper(embedding), "_1")) + ylab(paste0(toupper(embedding), "_2")) +
theme(plot.title = element_text(hjust = 0.5, margin=margin(0,0,0,0)),
panel.grid.minor = element_blank(), panel.grid.major = element_blank(),
legend.margin = margin(0,0,0,0), panel.background = element_blank())
# labels
if(label){
datax_group <- stats::aggregate(datax[,c("x","y")], list(datax[[group.by]]),mean)
colnames(datax_group) <- c(group.by, "x", "y")
g <- g +
geom_text(mapping = aes(x = .data[["x"]], y = .data[["y"]], label = .data[[group.by]]), datax_group, size = font.size)
}
# splitting
if(!is.null(split.by)){
if(is.null(nrow)){
unique_split <- unique(datax[[split.by]])
nrow <- ceiling(length(unique_split)/ncol)
}
g <- g +
facet_wrap(as.formula(paste("~", split.by)), nrow = nrow, ncol = ncol)
}
# return data
return(g)
}
####
## Embedding Feature Plot ####
####
#' vrEmbeddingFeaturePlot
#'
#' Plotting features of spatially resolved cells and spots on embeddings from multiple assays in a VoltRon object.
#'
#' @param object a VoltRon object
#' @param embedding the embedding type, i.e. pca, umap etc.
#' @param features a set of features to be visualized, either from \link{vrFeatures} of raw or normalized data or columns of the \link{Metadata}.
#' @param combine.features whether to combine all features in one plot
#' @param n.tile should points be aggregated into tiles before visualization (see \link{geom_tile}). Applicable only for cells and molecules
#' @param norm if TRUE, the normalized data is used
#' @param log if TRUE, data features (excluding metadata features) will be log transformed
#' @param assay assay name (exp: Assay1) or assay class (exp: Visium, Xenium), see \link{SampleMetadata}.
#' if NULL, the default assay will be used, see \link{vrMainAssay}.
#' @param ncol column wise number of plots, for \link{ggarrange}
#' @param nrow row wise number of plots, for \link{ggarrange}
#' @param font.size font size
#' @param pt.size point size
#' @param keep.scale whether unify all scales for all features or not
#' @param common.legend whether to use a common legend for all plots, see \link{ggarrange}
#' @param collapse.plots whether to combine all ggplots
#'
#' @import ggplot2
#' @importFrom dplyr arrange desc
#'
#' @export
vrEmbeddingFeaturePlot <- function(object, embedding = "pca", features = NULL, combine.features = FALSE, n.tile = 0, norm = TRUE, log = FALSE, assay = NULL, ncol = 2, nrow = NULL,
font.size = 2, pt.size = 1, keep.scale = "feature", common.legend = TRUE, collapse.plots = TRUE) {
# check object
if(!inherits(object, "VoltRon"))
stop("Please provide a VoltRon object!")
# features
if(is.null(features))
stop("You have to define at least one feature")
# sample metadata
sample.metadata <- SampleMetadata(object)
# get assay names
assay_names <- vrAssayNames(object, assay = assay)
# get entity type and metadata
metadata <- Metadata(object, assay = assay)
# get data
data_features <- features[features %in% vrFeatures(object, assay = assay)]
if(length(data_features) > 0){
normdata <- vrData(object, features = data_features, assay = assay_names, norm = norm)
if(log)
normdata <- log1p(normdata)
}
# get embedding
datax <- data.frame(vrEmbeddings(object, assay = assay_names, type = embedding))
datax <- datax[,seq_len(2)]
colnames(datax) <- c("x", "y")
# calculate limits for plotting, all for making one scale, feature for making multiple
limits <- Map(function(feat){
if(feat %in% vrFeatures(object, assay = assay)){
return(getRange(normdata[feat, ]))
} else {
if(feat %in% colnames(metadata)){
return(getRange(metadata[, feat]))
} else {
stop("Feature '", feat, "' cannot be found in data or metadata!")
}
}
}, features)
# configure titles
feature_title <- as.list(features)
names(feature_title) <- features
legend_title <- feature_title
# for each feature
i <- 1
gg <- list()
g <- ggplot()
all_data <- NULL
colors <- get_rasterization_colors(length(features))
for(feat in features){
# get data
if(feat %in% vrFeatures(object, assay = assay)){
if(inherits(normdata, "IterableMatrix")){
datax$score <- as.matrix(normdata[feat, rownames(datax)])[1,]
} else {
datax$score <- normdata[feat, rownames(datax)]
}
} else {
datax$score <- metadata[rownames(datax),feat]
}
# get image information and plotting features
midpoint <- sum(limits[[feat]])/2
# make one plot by combining features
if(combine.features){
# plot
g <- ggplot()
# add points, rasterize if requested or needed
if(n.tile > 0 || nrow(datax) > 50000){
if(n.tile == 0)
n.tile <- 1000
g <- g +
stat_summary_2d(mapping = aes(x = x, y = y, z = score), fun = mean, data = datax, geom = "tile", bins = n.tile, drop = TRUE) +
scale_fill_gradientn(name = legend_title[[feat]],
colors=c("grey97", colors[i]),
values=rescale_numeric(c(limits[[feat]][1], limits[[feat]][2])), limits = limits[[feat]])
all_data <- rbind(all_data,
data.frame(layer_data(g), color_group = colors[i]))
} else {
g <- g +
geom_point(mapping = aes(x = x, y = y, color = score), datax, shape = 16, size = pt.size) +
scale_color_gradientn(name = legend_title[[feat]],
colors=c("grey97", colors[i]),
values=rescale_numeric(c(limits[[feat]][1], limits[[feat]][2])), limits = limits[[feat]])
all_data <- rbind(all_data,
data.frame(layer_data(g), value = datax$score, color_group = colors[i]))
}
# make individual feature plots
} else {
# plot
g <- ggplot()
# add points, rasterize if requested or needed
if(n.tile > 0 || nrow(datax) > 50000){
if(n.tile == 0)
n.tile <- 1000
g <- vrFeaturePlotTiling(g = g, data = datax, legend_title = legend_title[[feat]], n.tile = n.tile, type = "embedding", limits = limits[[feat]])
} else {
g <- g +
geom_point(mapping = aes(x = x, y = y, color = score), dplyr::arrange(datax,score), shape = 16, size = pt.size) +
scale_color_gradientn(name = legend_title[[feat]],
colors=c("grey97", "blue"),
values=rescale_numeric(c(limits[[feat]][1], limits[[feat]][2])), limits = limits[[feat]])
}
}
# more visualization parameters
# g <- g + addEmbeddingPlotStyle(embedding)
g <- g +
theme_classic() +
xlab(paste0(toupper(embedding), "_1")) + ylab(paste0(toupper(embedding), "_2")) +
theme(plot.title = element_text(hjust = 0.5, margin=margin(0,0,0,0)),
panel.grid.minor = element_blank(), panel.grid.major = element_blank(),
legend.margin = margin(0,0,0,0), panel.background = element_blank())
# add graph to list
gg[[i]] <- g
i <- i + 1
}
# collapse plots
if(collapse.plots){
# return a list of plots or a single one
if(length(features) > 1){
# combine features
if(combine.features){
g.combined <- vrEmbeddingFeatureCombinePlot(all_data, n.tile, datax, features, embedding)
return(g.combined)
} else {
if(length(gg) < ncol) ncol <- length(gg)
return(ggpubr::ggarrange(plotlist = gg, ncol = ncol, nrow = ceiling(length(gg)/ncol)))
}
} else {
return(gg[[1]])
}
} else {
return(gg)
}
}
#' vrEmbeddingFeatureCombinePlot
#'
#' A different plotting setting for vrEmbeddingFeaturePlot when combine.feature=TRUE
#'
#' @param all_data summary data
#' @param n.tile should points be aggregated into tiles before visualization (see \link{geom_tile}). Applicable only for cells and molecules
#' @param datax original plotting data
#' @param features features
#' @param embedding embedding
#'
#' @import ggplot2
#'
#' @noRd
vrEmbeddingFeatureCombinePlot <- function(all_data, n.tile, datax, features, embedding){
# ggplot
g.combined <- ggplot()
# tiling or not
if(n.tile > 0 || nrow(datax) > 50000){
if(n.tile == 0)
n.tile <- 1000
all_data <- all_data %>% group_by(x,y) %>%
summarize(fill = fill[which.max(value)],
group = color_group[which.max(value)],
value = value[which.max(value)])
key_table <- all_data[,c("fill", "group", "value")] %>%
dplyr::group_by(group) %>%
dplyr::summarise(fill = fill[which.max(value)], value = max(value))
g.combined <- g.combined +
ggplot2::geom_tile(data = as.data.frame(all_data), aes(x = x, y = y, fill = fill)) +
ggplot2::scale_fill_identity("", labels = features, breaks = key_table$fill, guide = "legend")
} else {
all_data <- all_data %>% group_by(x,y) %>%
dplyr::summarise(color = colour[which.max(value)],
group = color_group[which.max(value)],
value = value[which.max(value)])
key_table <- all_data[,c("color", "group", "value")] %>%
dplyr::group_by(group) %>%
dplyr::summarise(color = color[which.max(value)], value = max(value))
g.combined <- g.combined +
ggplot2::geom_point(data = as.data.frame(all_data), aes(x = x, y = y, color = color)) +
ggplot2::scale_color_identity("", labels = features, breaks = key_table$color, guide = "legend")
}
# style
g.combined <- g.combined +
addEmbeddingPlotStyle(embedding) +
xlab(paste0(toupper(embedding), "_1")) +
ylab(paste0(toupper(embedding), "_2"))
g.combined
}
####
# Scatter Plot ####
####
#' vrScatterPlot
#'
#' get a scatter plot between two features
#'
#' @param object a VoltRon object
#' @param feature.1 first feature
#' @param feature.2 second feature
#' @param norm if TRUE, the normalize data will be used
#' @param assay assay name (exp: Assay1) or assay class (exp: Visium, Xenium), see \link{SampleMetadata}.
#' if NULL, the default assay will be used, see \link{vrMainAssay}.
#' @param pt.size point size
#' @param font.size font size
#' @param group.by a column of metadata from \link{Metadata} used as grouping label for the spatial entities
#' @param label whether labels are visualized or not
#' @param trend inserting a trend line two the scatter plot
#'
#' @import ggplot2
#' @importFrom ggrepel geom_label_repel
#'
#' @export
#'
vrScatterPlot <- function(object, feature.1, feature.2, norm = TRUE, assay = NULL,
pt.size = 2, font.size = 2, group.by = "label", label = FALSE, trend = FALSE){
# check object
if(!inherits(object, "VoltRon"))
stop("Please provide a VoltRon object!")
# check the number of features
if(is.null(feature.1) | is.null(feature.2))
stop("Please provide both 'feature.1' and 'feature.2'")
# check the number of features
if((length(feature.1) != 1 | length(feature.2) != 1))
stop("Both 'feature.1' and 'feature.2' should be of length 1.")
# data
normdata <- vrData(object, assay = assay, norm = norm)
# sample metadata
sample.metadata <- SampleMetadata(object)
# get assay names
assay_names <- vrAssayNames(object, assay = assay)
# get entity type and metadata
metadata <- Metadata(object, assay = assay)
# get data
# data_feature <- sapply(c(feature.1, feature.2), function(feat){
data_feature <- vapply(c(feature.1, feature.2), function(feat){
if(feat %in% rownames(normdata)){
return(normdata[feat,])
} else {
return(metadata[,feat])
}
}, numeric(nrow(metadata)))
data_feature <- as.data.frame(data_feature)
# plot
g <- ggplot()
# plot scatter
g <- g +
geom_point(mapping = aes(x = .data[[feature.1]], y = .data[[feature.2]]), data = data_feature, size = pt.size)
# visualize labels
if(label){
data_feature[[group.by]] <- metadata[,group.by]
g <- g + ggrepel::geom_label_repel(mapping = aes(x = .data[[feature.1]], y = .data[[feature.2]], label = .data[[group.by]]), data_feature,
box.padding = 0.5, size = font.size, direction = "both", seed = 1)
}
# visualize trend
if(trend){
g <- g + geom_smooth()
}
g
}
####
# Heatmap Plot ####
####
#' HeatmapPlot
#'
#' @param object a VoltRon object
#' @param assay assay name (exp: Assay1) or assay class (exp: Visium, Xenium), see \link{SampleMetadata}.
#' if NULL, the default assay will be used, see \link{vrMainAssay}.
#' @param features a set of features to be visualized from \link{vrFeatures} of raw or normalized data.
#' @param group.by a column of metadata from \link{Metadata} used as grouping label for the spatial entities
#' @param norm if TRUE, the normalized data is used
#' @param scaled if TRUE, the data will be scaled before visualization
#' @param show_row_names if TRUE, row names of the heatmap will be shown
#' @param cluster_rows if TRUE, the rows of the heatmap will be clustered
#' @param show_heatmap_legend if TRUE, the heatmap legend is shown
#' @param outlier.quantile quantile for detecting outliers whose values are set to the quantile, change to lower values to adjust large number of outliers, default: 0.99
#' @param highlight.some if TRUE, some rows will be showed at random, reproducible by \code{seed} argument
#' @param n_highlight the number of row labels shown, if \code{show_row_names} is TRUE
#' @param font.size font size
#' @param ... additional parameters passed to \link{getVariableFeatures}
#'
#' @importFrom stats quantile
#'
#' @export
vrHeatmapPlot <- function(object, assay = NULL, features = NULL, group.by = "clusters",
norm = TRUE, scaled = TRUE, show_row_names = NULL, cluster_rows = TRUE, show_heatmap_legend = FALSE,
outlier.quantile = 0.99, highlight.some = FALSE, n_highlight = 30, font.size = 13.2, ...){
if (!requireNamespace('ComplexHeatmap'))
stop("Please install ComplexHeatmap package to use the Heatmap function!: BiocManager::install('ComplexHeatmap')")
if (!requireNamespace('viridisLite'))
stop("Please install viridisLite package to use the Heatmap function!: install.packages('viridisLite')")
# check object
if(!inherits(object, "VoltRon"))
stop("Please provide a VoltRon object!")
# data
heatmapdata <- vrData(object, assay = assay, norm = norm)
# features
if(is.null(features)){
if(nrow(vrFeatureData(object, assay = assay)) > 0){
features <- getVariableFeatures(object, assay = assay, ...)
} else {
features <- vrFeatures(object, assay = assay)
}
} else {
nonmatching_features <- setdiff(features, vrFeatures(object, assay = assay))
features <- features[features %in% vrFeatures(object, assay = assay)]
if(length(features) == 0){
stop("None of the provided features are found in the assay")
}
if(length(nonmatching_features))
message("the following features are not found in the assay: ", paste(nonmatching_features, collapse = ", "))
}
heatmapdata <- heatmapdata[features, ]
# get entity type and metadata
metadata <- Metadata(object, assay = assay)
if(!is.null(rownames(metadata))){
metadata <- metadata[colnames(heatmapdata),]
} else {
metadata <- metadata[match(colnames(heatmapdata), as.vector(metadata$id)),]
}
# scaling, optional
if(scaled){
heatmapdata <- apply(heatmapdata, 1, scale)
heatmapdata <- t(heatmapdata)
legend_title <- "Scaled \n Exp."
} else {
legend_title <- "Norm. \n Exp."
}
# check for nan
heatmapdata[is.nan(heatmapdata)] <- 0
# manage data for plotting
if(group.by %in% colnames(metadata)){
heatmapdata <- heatmapdata[,order(metadata[[group.by]], decreasing = FALSE)]
labels_ordered <- metadata[[group.by]][order(metadata[[group.by]], decreasing = FALSE)]
labels_ordered_table <- table(labels_ordered)
col_split = factor(labels_ordered, levels = names(labels_ordered_table))
} else {
stop("Column '", group.by, "' cannot be found in metadata!")
}
# update limits
limits <- stats::quantile(as.vector(heatmapdata), probs = c(1-outlier.quantile, outlier.quantile))
heatmapdata[heatmapdata > limits[2]] <- limits[2]
heatmapdata[heatmapdata < limits[1]] <- limits[1]
# highlight some rows
if(highlight.some){
ind <- sample(seq_len(nrow(heatmapdata)), n_highlight, replace = FALSE)
ha <- ComplexHeatmap::rowAnnotation(foo = ComplexHeatmap::anno_mark(at = ind,
labels = rownames(heatmapdata)[ind],
padding = 1,
labels_gp = gpar(fontsize = font.size)))
} else{
ha <- NULL
}
# visualize
if(is.null(show_row_names))
show_row_names <- (nrow(heatmapdata) < 30)
legend_at <- seq(min(heatmapdata), max(heatmapdata), (max(heatmapdata)-min(heatmapdata))/5)
legend_label <- round(legend_at, 2)
ComplexHeatmap::Heatmap(heatmapdata,
show_row_names = show_row_names, show_row_dend = FALSE, row_names_gp = gpar(fontsize = font.size),
show_column_names = FALSE, column_title_rot = 45, column_title_gp = gpar(fontsize = font.size),
column_split = col_split, cluster_columns = FALSE, cluster_rows = cluster_rows,
show_heatmap_legend = show_heatmap_legend,
heatmap_legend_param = list(title = legend_title, at = legend_at, labels = legend_label),
right_annotation = ha,
col = viridisLite::viridis(100))
}
####
# Violin Plot ####
####
#' vrViolinPlot
#'
#' @param object a VoltRon object
#' @param features a set of features to be visualized, either from \link{vrFeatures} of raw or normalized data or columns of the \link{Metadata}.
#' @param assay assay name (exp: Assay1) or assay class (exp: Visium, Xenium), see \link{SampleMetadata}.
#' if NULL, the default assay will be used, see \link{vrMainAssay}.
#' @param group.by a column of metadata from \link{Metadata} used as grouping label for the spatial entities
#' @param norm if TRUE, the normalized data is used
#' @param pt.size point size
#' @param plot.points if TRUE, measures are visualized as points as well.
#' @param ncol column wise number of plots, for \link{ggarrange}
#' @param nrow row wise number of plots, for \link{ggarrange}
#'
#' @import ggplot2
#' @importFrom data.table data.table melt
#'
#' @export
#'
vrViolinPlot <- function(object, features = NULL, assay = NULL, group.by = "Sample",
norm = TRUE, pt.size = 0.5, plot.points = TRUE, ncol = 2, nrow = NULL){
# check object
if(!inherits(object, "VoltRon"))
stop("Please provide a VoltRon object!")
# features
if(is.null(features))
stop("You have to define at least one feature")
# sample metadata
sample.metadata <- SampleMetadata(object)
# get entity type and metadata
metadata <- Metadata(object, assay = assay)
# get data
if(any(features %in% vrFeatures(object, assay = assay))){
selected_features <- features[features %in% vrFeatures(object, assay = assay)]
violindata <- vrData(object, features = selected_features, assay = assay, norm = norm)
}
# get feature data
datax <- lapply(features, function(x){
if(x %in% vrFeatures(object, assay = assay)){
return(violindata[x,])
} else if(x %in% colnames(metadata)){
return(metadata[,x])
} else {
stop("Please provide feature names which reside in either the data or metadata slots!")
}
})
datax <- do.call(cbind, datax)
colnames(datax) <- features
# assays
assays <- stringr::str_extract(rownames(metadata), "Assay[0-9]+$")
assay_title <- apply(sample.metadata[assays,], 1, function(x) paste(x["Sample"], x["Layer"], x["Assay"], sep = "|"))
# violin plot
ggplotdatax <- data.frame(datax,
group.by = as.factor(metadata[[group.by]]),
assay_title = assay_title,
spatialpoints = rownames(metadata))
# ggplotdatax <- reshape2::melt(ggplotdatax, id.var = c("group.by", "assay_title", "spatialpoints"))
ggplotdatax <- data.table::melt(data.table::data.table(ggplotdatax), id.var = c("group.by", "assay_title", "spatialpoints"))
# visualize points on violin
if(plot.points){
gg <- ggplot(ggplotdatax, aes(x = group.by, y = value, color = group.by)) +
geom_violin() +
geom_point(size = pt.size, position = position_jitter()) +
guides(fill = guide_legend(show = FALSE), color = guide_legend(title = group.by, override.aes=list(size = 2)))
} else {
gg <- ggplot(ggplotdatax, aes(x = group.by, y = value, color = group.by, fill = group.by)) +
geom_violin() +
guides(color = guide_legend(title = group.by, override.aes=list(size = 2)), fill = guide_legend(title = group.by, override.aes=list(size = 2)))
}
# theme
gg <- gg +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
ylab("") + xlab(group.by)
if(length(features) > 1){
if(length(gg) < ncol) ncol <- length(gg)
gg <- gg + facet_wrap(.~variable, ncol = ncol, nrow = ceiling(length(gg)/ncol), scales = "free_y")
} else {
gg <- gg + labs(title = features)
}
return(gg)
}
####
# ROI Plots ####
####
#' vrBarPlot
#'
#' @param object a VoltRon object
#' @param features a set of features to be visualized, either from \link{vrFeatures} of raw or normalized data or columns of the \link{Metadata}.
#' @param assay assay name (exp: Assay1) or assay class (exp: Visium, Xenium), see \link{SampleMetadata}.
#' if NULL, the default assay will be used, see \link{vrMainAssay}.
#' @param x.label labels of the x axis
#' @param group.by a column of metadata from \link{Metadata} used as grouping label for the spatial entities
#' @param split.by the column to split the barplots by
#' @param norm if TRUE, the normalized data is used
#' @param log if TRUE, data features (excluding metadata features) will be log transformed
#' @param ncol column wise number of plots, for \link{ggarrange}
#' @param nrow row wise number of plots, for \link{ggarrange}
#'
#' @import ggplot2
#' @importFrom data.table data.table melt
#'
#' @export
vrBarPlot <- function(object, features = NULL, assay = NULL, x.label = NULL, group.by = "Sample",
split.by = NULL, norm = TRUE, log = FALSE, ncol = 2, nrow = NULL){
# check object
if(!inherits(object, "VoltRon"))
stop("Please provide a VoltRon object!")
# features
if(is.null(features))
stop("You have to define at least one feature")
# sample metadata
sample.metadata <- SampleMetadata(object)
# get assay names
assay_names <- vrAssayNames(object, assay = assay)
# get data
barplotdata <- vrData(object, assay = assay, norm = norm)
if(log)
barplotdata <- log1p(barplotdata)
# get entity type and metadata
assay_types <- vrAssayTypes(object, assay = assay)
if(unique(assay_types) %in% c("spot","cell")){
stop("vrBarPlot can only be used for ROI assays")
} else {
metadata <- Metadata(object, type = "ROI")
metadata <- subset_metadata(metadata, assays = assay_names)
}
# get feature data
datax <- lapply(features, function(x){
if(x %in% rownames(barplotdata)){
return(as.vector(as(barplotdata[x,,drop = FALSE], "dgCMatrix")))
} else if(x %in% colnames(metadata)){
return(as.vector(metadata[,x]))
} else{
stop("Feature '", x, "' cannot be found in data or metadata!")
}
})
datax <- as.data.frame(do.call(cbind, datax))
colnames(datax) <- features
# get assay titles
assays <- stringr::str_extract(vrSpatialPoints(object, assay = assay_names), "Assay[0-9]+$")
assay_title <- apply(sample.metadata[assays,], 1, function(x) paste(x["Sample"], x["Layer"], x["Assay"], sep = "|"))
# labels and groups
if(is.null(x.label)) {
if(is.null(rownames(metadata))){
x.labels <- factor(metadata$id)
} else {
x.labels <- factor(rownames(metadata))
}
} else {
if(x.label %in% colnames(metadata)){
x.labels <- factor(as.vector(metadata[[x.label]]))
} else {
stop("Column '", x.label, "' cannot be found in metadata!")
}
}
if(group.by %in% colnames(metadata)){
group.by.col <- factor(as.vector(metadata[[group.by]]))
} else {
stop("Column '", group.by, "' cannot be found in metadata!")
}
# plotting data
if(!is.null(rownames(metadata))){
spatialpoints <- rownames(metadata)
} else {
spatialpoints <- metadata$id
}
if(is.null(split.by)){
ggplotdatax <- data.frame(datax,
x.label = x.labels,
group.by = group.by.col,
assay_title = assay_title,
spatialpoints = spatialpoints)
ggplotdatax <- data.table::melt(data.table::data.table(ggplotdatax), id.var = c("x.label", "assay_title", "group.by", "spatialpoints"))
gg <- ggplot(ggplotdatax, aes(x = x.label, y = value,
fill = factor(group.by, levels = unique(group.by)))) +
geom_bar(stat = "identity") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
ylab("") + xlab(x.label) +
guides(fill = guide_legend(title = group.by))
if(length(features) > 1){
if(length(gg) < ncol) ncol <- length(gg)
gg <- gg + facet_wrap(.~variable, ncol = ncol, nrow = ceiling(length(features)/ncol), scales = "free")
return(gg)
} else {
gg <- gg + labs(title = features)
return(gg)
}
} else {
# check split.by column name
if(split.by %in% colnames(metadata)){
split.by.col <- factor(as.vector(metadata[[split.by]]))
} else {
stop("Column '", split.by, "' cannot be found in metadata!")
}
# make ggplot
ggplotdatax <- data.frame(datax,
x.label = x.labels,
group.by = group.by.col,
split.by = split.by.col,
assay_title = assay_title,
spatialpoints = spatialpoints)
ggplotdatax <- data.table::melt(data.table::data.table(ggplotdatax), id.var = c("x.label", "assay_title", "group.by", "split.by", "spatialpoints"))
gg <- ggplot(ggplotdatax, aes(x = x.label, y = value,
fill = factor(group.by, levels = unique(group.by)))) +
geom_bar(stat = "identity") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
ylab("") + xlab(x.label) +
guides(fill = guide_legend(title = group.by))
if(length(features) > 1){
gg <- gg + facet_grid(variable~split.by, scales = "free", space = "free_x")
return(gg)
} else {
if(length(gg) < ncol) ncol <- length(gg)
gg <- gg + facet_wrap(.~split.by, ncol = ncol, nrow = ceiling(length(unique(split.by.col))/ncol), scales = "free_x")
return(gg)
}
}
}
#' vrPercentagePlot
#'
#' @param object a VoltRon object
#' @param assay assay name (exp: Assay1) or assay class (exp: Visium, Xenium), see \link{SampleMetadata}.
#' if NULL, the default assay will be used, see \link{vrMainAssay}.
#' @param x.label labels of the x axis
#' @param split.by the column to split the barplots by
#' @param split.method either facet_grid or facet_wrap, not affected if \code{split.by} is \code{NULL}
#' @param ncol column wise number of plots, for \link{ggarrange}
#' @param nrow row wise number of plots, for \link{ggarrange}
#'
#' @import ggplot2
#' @importFrom data.table data.table melt
#'
#' @export
#'
vrProportionPlot <- function(object, assay = NULL, x.label = NULL,
split.by = NULL, split.method = "facet_wrap", ncol = 2, nrow = NULL){
# check object
if(!inherits(object, "VoltRon"))
stop("Please provide a VoltRon object!")
# sample metadata
sample.metadata <- SampleMetadata(object)
# get assay names
assay_names <- vrAssayNames(object, assay = assay)
# get data
barplotdata <- vrData(object, assay = assay, norm = FALSE)
# get entity type and metadata
assay_types <- vrAssayTypes(object, assay = assay)
if(unique(assay_types) %in% c("spot","cell")){
stop("vrProportionPlot can only be used for ROI assays")
} else {
metadata <- Metadata(object, type = "ROI")
metadata <- subset_metadata(metadata, assays = assay_names)
}
# violin plot
assays <- stringr::str_extract(vrSpatialPoints(object, assay = assay_names), "Assay[0-9]+$")
assay_title <- apply(sample.metadata[assays,], 1, function(x) paste(x["Sample"], x["Layer"], x["Assay"], sep = "|"))
# labels and groups
if(is.null(x.label)) {
if(is.null(rownames(metadata))){
x.labels <- factor(metadata$id)
} else {
x.labels <- factor(rownames(metadata))
}
} else {
if(x.label %in% colnames(metadata)){
x.labels <- factor(as.vector(metadata[[x.label]]))
} else {
stop("Column '", x.label, "' cannot be found in metadata!")
}
}
# plotting data
if(!is.null(rownames(metadata))){
spatialpoints <- rownames(metadata)
} else {
spatialpoints <- metadata$id
}
ggplotdatax <- data.frame(t(barplotdata),
x.label = x.label,
assay_title = assay_title,
spatialpoints = spatialpoints)
ggplotdatax <- data.table::melt(data.table::data.table(ggplotdatax), id.var = c("x.label", "assay_title", "spatialpoints"))
ggplotdatax <- ggplotdatax[ggplotdatax$value > 0,]
gg <- ggplot(ggplotdatax, aes(x = x.label, y = value, fill = variable)) +
geom_bar(stat = "identity") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
ylab("") + xlab("") +
guides(fill = guide_legend(title = ""))
if(is.null(split.by)){
ggplotdatax <- data.frame(t(barplotdata),
x.label = x.label,
assay_title = assay_title,
spatialpoints = rownames(metadata))
ggplotdatax <- data.table::melt(data.table::data.table(ggplotdatax), id.var = c("x.label", "assay_title", "spatialpoints"))
ggplotdatax <- ggplotdatax[ggplotdatax$value > 0,]
gg <- ggplot(ggplotdatax, aes(x = x.label, y = value, fill = variable)) +
geom_bar(stat = "identity") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
ylab("") + xlab("") +
guides(fill = guide_legend(title = ""))
} else {
# check split.by column name
if(split.by %in% colnames(metadata)){
split.by.col <- factor(metadata[[split.by]])
} else {
stop("Column '", split.by, "' cannot be found in metadata!")
}
ggplotdatax <- data.frame(t(barplotdata),
x.label = x.label,
assay_title = assay_title,
split.by = split.by.col,
spatialpoints = spatialpoints)
ggplotdatax <- data.table::melt(data.table::data.table(ggplotdatax), id.var = c("x.label", "assay_title", "split.by", "spatialpoints"))
ggplotdatax <- ggplotdatax[ggplotdatax$value > 0,]
gg <- ggplot(ggplotdatax, aes(x = x.label, y = value, fill = variable)) +
geom_bar(stat = "identity") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
ylab("") + xlab("") +
guides(fill = guide_legend(title = ""))
# split
if(length(gg) < ncol) ncol <- length(gg)
if(split.method == "facet_wrap"){
gg <- gg + facet_wrap(.~split.by, ncol = ncol, nrow = ceiling(length(unique(split.by.col))/ncol), scales = "free_x")
} else {
gg <- gg + facet_grid(.~split.by, scales = "free_x", space = "free_x", )
}
}
return(gg)
}
####
# Rasterization ####
####
#' vrGroupPlotTiling
#'
#' Plotting a tiled version of the vrEmbeddingPlot and vrSpatialPlot
#'
#' @param g the ggplot figure
#' @param data the data frame with coordinates and group identities
#' @param group.by a column of metadata from \link{Metadata} used as grouping label for the spatial entities
#' @param n.tile should points be aggregated into tiles before visualization (see \link{geom_tile}). Applicable only for cells and molecules
#' @param alpha alpha level of colors of visualized points and segments
#'
#' @import ggplot2
#'
#' @noRd
vrGroupPlotTiling <- function(g, data, group.by, n.tile, alpha = 1){
g + stat_bin_2d(mapping = aes(x = .data[["x"]], y = .data[["y"]], fill = .data[[group.by]], color = .data[[group.by]]),
data = data, bins = n.tile, drop = TRUE, alpha = alpha)
}
#' vrFeaturePlotTiling
#'
#' Plotting a tiled version of the vrEmbeddingFeaturePlot and vrSpatialFeaturePlot
#'
#' @param g the ggplot figure
#' @param data the data frame with coordinates and group identities
#' @param legend_title the legend title of the single plot
#' @param n.tile should points be aggregated into tiles before visualization (see \link{geom_tile}). Applicable only for cells and molecules
#' @param alpha alpha level of colors of visualized points and segments
#' @param type embedding or spatial, changes the color scheme
#'
#' @import ggplot2
#'
#' @noRd
vrFeaturePlotTiling <- function(g, data, legend_title, n.tile, alpha = 1, limits = NULL, type = "embedding"){
# get summary per title
gplot <- g + stat_summary_2d(mapping = aes(x = x, y = y, z = score), fun = mean,
data = data, geom = "tile", bins = n.tile, drop = TRUE, alpha = alpha)
# check if limits are defined, typically coming from vrEmbeddingFeaturePlot
if(is.null(limits)){
hex_count_data <- suppressWarnings(ggplot_build(gplot)$data)
hex_count_data <- hex_count_data[[length(hex_count_data)]]
midpoint <- max(hex_count_data$value)/2
}
# color scheme for either spatial or embedding feature plot
if(type == "spatial"){
gplot <- gplot +
scale_fill_gradientn(name = legend_title,
colors=c("dodgerblue2", "white", "yellow3"),
values=rescale_numeric(c(min(hex_count_data$value), midpoint, max(hex_count_data$value)),
limits = c(min(hex_count_data$value), max(hex_count_data$value))))
} else{
gplot <- gplot +
scale_fill_gradientn(name = legend_title,
colors=c("grey97", "blue"),
values=rescale_numeric(c(limits[1], limits[2])), limits = limits)
}
# return
gplot
}
#' @param n number of colors
#'
#' @noRd
get_rasterization_colors <- function(n){
colors <- c("darkblue",'darkred', "yellow2", "darkgreen", "darkorange", 'darkmagenta', "brown")
colors[seq_len(n)]
}
####
# Plotting style
####
# this function is deprecated
addEmbeddingPlotStyle <- function(embedding){
suppressWarnings(
theme_classic() +
xlab(paste0(toupper(embedding), "_1")) + ylab(paste0(toupper(embedding), "_2")) +
theme(plot.title = element_text(hjust = 0.5, margin=margin(0,0,0,0)),
panel.grid.minor = element_blank(), panel.grid.major = element_blank(),
legend.margin = margin(0,0,0,0), panel.background = element_blank())
)
}
Datasets
Models
