diff --git a/DESCRIPTION b/DESCRIPTION index a650e263..85b69604 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: miaViz Title: Microbiome Analysis Plotting and Visualization -Version: 1.21.1 +Version: 1.21.2 Authors@R: c(person(given = "Tuomas", family = "Borman", role = c("aut", "cre"), email = "tuomas.v.borman@utu.fi", @@ -31,7 +31,7 @@ License: Artistic-2.0 | file LICENSE Encoding: UTF-8 LazyData: false Depends: - R (>= 4.0), + R (>= 4.1), ggplot2, ggraph (>= 2.0), mia (>= 1.13.0), @@ -47,6 +47,7 @@ Imports: ggnewscale, ggrepel, ggtree, + grid, methods, patchwork, rlang, diff --git a/NAMESPACE b/NAMESPACE index 68b4e23b..fa0a1614 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -10,9 +10,11 @@ export(plotDMNFit) export(plotFeaturePrevalence) export(plotForest) export(plotHistogram) +export(plotJointRPCA) export(plotLoadings) export(plotMediation) export(plotNMDS) +export(plotOrdination) export(plotPrevalence) export(plotPrevalentAbundance) export(plotRDA) @@ -39,8 +41,10 @@ exportMethods(plotDMNFit) exportMethods(plotFeaturePrevalence) exportMethods(plotForest) exportMethods(plotHistogram) +exportMethods(plotJointRPCA) exportMethods(plotLoadings) exportMethods(plotMediation) +exportMethods(plotOrdination) exportMethods(plotPrevalence) exportMethods(plotPrevalentAbundance) exportMethods(plotRDA) @@ -68,7 +72,6 @@ importFrom(BiocParallel,bpstop) importFrom(DelayedArray,rowMeans) importFrom(DelayedArray,rowSums) importFrom(DirichletMultinomial,mixture) -importFrom(S4Vectors,SimpleList) importFrom(S4Vectors,metadata) importFrom(S4Vectors,unfactor) importFrom(SingleCellExperiment,reducedDim) @@ -80,6 +83,7 @@ importFrom(SummarizedExperiment,colData) importFrom(SummarizedExperiment,rowData) importFrom(TreeSummarizedExperiment,colTree) importFrom(TreeSummarizedExperiment,colTreeNames) +importFrom(TreeSummarizedExperiment,rowLinks) importFrom(TreeSummarizedExperiment,rowTree) importFrom(TreeSummarizedExperiment,rowTreeNames) importFrom(ape,as.phylo) @@ -106,6 +110,7 @@ importFrom(dplyr,rename) importFrom(dplyr,rename_with) importFrom(dplyr,row_number) importFrom(dplyr,select) +importFrom(dplyr,slice_max) importFrom(dplyr,summarise) importFrom(dplyr,summarize) importFrom(dplyr,ungroup) @@ -147,8 +152,6 @@ importFrom(ggplot2,scale_y_discrete) importFrom(ggplot2,theme) importFrom(ggplot2,theme_bw) importFrom(ggplot2,theme_classic) -importFrom(ggplot2,xlab) -importFrom(ggplot2,ylab) importFrom(ggrepel,geom_label_repel) importFrom(ggrepel,geom_text_repel) importFrom(ggtree,geom_cladelab) @@ -162,6 +165,8 @@ importFrom(ggtree,ggtree) importFrom(ggtree,gheatmap) importFrom(ggtree,groupOTU) importFrom(ggtree,theme_tree) +importFrom(grid,arrow) +importFrom(grid,unit) importFrom(mia,meltSE) importFrom(patchwork,plot_layout) importFrom(patchwork,wrap_plots) @@ -170,7 +175,6 @@ importFrom(rlang,":=") importFrom(rlang,sym) importFrom(scales,pretty_breaks) importFrom(scales,rescale) -importFrom(scater,plotReducedDim) importFrom(scater,retrieveCellInfo) importFrom(scater,retrieveFeatureInfo) importFrom(stats,formula) diff --git a/NEWS b/NEWS index b55174bf..aa5b7d8a 100644 --- a/NEWS +++ b/NEWS @@ -51,3 +51,6 @@ Changes in version 1.17.x Changes in version 1.19.x + plotRDA: Now plotting works with interaction term (2025-11-09) + plotBoxplot: Added option to add p-values (2026-01-07) + +Changes in version 1.21.x ++ Added plotOrdination and plotJointRPCA (2026-08-04) diff --git a/R/AllGenerics.R b/R/AllGenerics.R index 527c39a0..1058ba59 100644 --- a/R/AllGenerics.R +++ b/R/AllGenerics.R @@ -156,3 +156,13 @@ setGeneric("plotBarplot", signature = c("x"), function(x, ...) #' @export setGeneric("plotBoxplot", signature = c("object"), function(object, ...) standardGeneric("plotBoxplot")) + +#' @rdname plotOrdination +#' @export +setGeneric("plotOrdination", signature = c("x"), function(x, ...) + standardGeneric("plotOrdination")) + +#' @rdname plotJointRPCA +#' @export +setGeneric("plotJointRPCA", signature = c("x"), function(x, ...) + standardGeneric("plotJointRPCA")) diff --git a/R/plotBoxplot.R b/R/plotBoxplot.R index 4e7d739d..89ab1739 100644 --- a/R/plotBoxplot.R +++ b/R/plotBoxplot.R @@ -1106,14 +1106,16 @@ setMethod("plotBoxplot", signature = c(object = "SummarizedExperiment"), # This function adds points to plot .add_points_layer <- function( - p, df, point.alpha = 0.65, point.size = 2, point.shape = 19L, + p, df, + x = NULL, y = NULL, + point.alpha = 0.65, point.size = 2, point.shape = 19L, point.colour = point.color, point.color = "grey70", ...){ # To disable "no visible binding for global variable" message in cmdcheck x_point <- y_point <- NULL args <- list( mapping = aes( - x = x_point, - y = y_point, + x = if(!is.null(x)) .data[[x]] else x_point, + y = if(!is.null(y)) .data[[y]] else y_point, colour = if(!is.null(attributes(df)[["colour.by"]])) .data[[attributes(df)[["colour.by"]]]], shape = if(!is.null(attributes(df)[["shape.by"]])) @@ -1332,5 +1334,22 @@ setMethod("plotBoxplot", signature = c(object = "SummarizedExperiment"), if( !is.null(attributes(df)[["size.by"]]) ){ p <- p + labs(shape = attributes(df)[["size.by"]]) } + + # Adjust colors + if( !is.null(attributes(df)[["colour.by"]]) ){ + name <- attributes(df)[["colour.by"]] + vals <- df[[name]] + p <- .resolve_plot_colours( + p, vals, name, + ) + } + if( !is.null(attributes(df)[["fill.by"]]) ){ + name <- attributes(df)[["fill.by"]] + vals <- df[[name]] + p <- .resolve_plot_colours( + p, vals, name, fill = TRUE + ) + } + return(p) } diff --git a/R/plotCCA.R b/R/plotCCA.R index 0bc0fbb7..80c851ed 100644 --- a/R/plotCCA.R +++ b/R/plotCCA.R @@ -6,8 +6,7 @@ #' #' @param x a #' \code{\link[TreeSummarizedExperiment:TreeSummarizedExperiment-constructor]{TreeSummarizedExperiment}} -#' or a matrix of weights. The latter is returned as output from -#' \code{\link[mia:runCCA]{getRDA}}. +#' object. #' #' @param dimred \code{Character scalar} or \code{integer scalar}. Determines #' the reduced dimension to @@ -119,14 +118,17 @@ #' formula = assay ~ ClinicalStatus + Gender + Age, #' distance = "bray", #' na.action = na.exclude -#' ) +#' ) #' #' suppressWarnings({ #' # Create RDA plot coloured by variable #' plotRDA(tse, "RDA", colour.by = "ClinicalStatus") #' -#' # Create RDA plot with empty ellipses -#' plotRDA(tse, "RDA", colour.by = "ClinicalStatus", add.ellipse = "colour") +#' # Create RDA plot with ellipses +#' plotRDA( +#' tse, "RDA", colour.by = "ClinicalStatus", fill.by = "ClinicalStatus", +#' add.ellipse = TRUE +#' ) #' #' # Create RDA plot with text encased in labels #' plotRDA(tse, "RDA", colour.by = "ClinicalStatus", vec.text = FALSE) @@ -136,20 +138,13 @@ #' #' # Create RDA plot without vectors #' plotRDA(tse, "RDA", colour.by = "ClinicalStatus", add.vectors = FALSE) -#' -#' # Calculate RDA as a separate object -#' rda_mat <- getRDA( -#' tse, -#' assay.type = "relabundance", -#' formula = assay ~ ClinicalStatus + Gender + Age, -#' distance = "bray", -#' na.action = na.exclude -#' ) -#' -#' # Create RDA plot from RDA matrix -#' plotRDA(rda_mat) #' }) #' +#' @seealso +#' \itemize{ +#' \item \code{\link[=plotOrdination]{plotOrdination}} +#' } +#' NULL #' @rdname plotCCA @@ -162,16 +157,6 @@ setMethod("plotCCA", signature = c(x = "SingleCellExperiment"), } ) -#' @rdname plotCCA -#' @aliases plotRDA -#' @export -setMethod("plotCCA", signature = c(x = "matrix"), - function(x, ...){ - # Reproduce plotRDA function - return(plotRDA(x, ...)) - } -) - #' @rdname plotCCA #' @aliases plotCCA #' @export @@ -198,73 +183,16 @@ setMethod("plotRDA", signature = c(x = "SingleCellExperiment"), args[["ncomponents"]] <- 2L # Get data for plotting plot_args <- list() - plot_args[["ellipse_data"]] <- do.call(.get_rda_ellipse_data, args) plot_args[["vector_data"]] <- do.call(.get_rda_vector_data, args) - plot_args[["centroids"]] <- do.call(.get_rda_centroids_data, args) - plot_args[["species_scores"]] <- do.call(.get_rda_species_data, args) - plot_args[["plot"]] <- do.call(.create_rda_baseplot, args) + p <- plotOrdination(x, dimred = dimred, ...) # Create a final plot - p <- .rda_plotter(plot_args, ...) - return(p) - } -) - -#' @rdname plotCCA -#' @aliases plotCCA -#' @export -setMethod("plotRDA", signature = c(x = "matrix"), - function(x, ...){ - # Construct TreeSE from rda/cca object - x <- .rda2tse(x) - # Run plotRDA method for TreeSE - p <- plotRDA(x, "RDA", ...) + p <- .rda_plotter_vector(p, plot_args, ...) return(p) } ) ################################ HELP FUNCTIONS ################################ -# Construct TreeSE from matrix to pass it to downstream functions. It is useful -# for instance if get* functios was used instead of add*. -#' @importFrom S4Vectors SimpleList -.rda2tse <- function(object) { - # Convert rda/cca object to TreeSE - object <- TreeSummarizedExperiment( - assays = SimpleList(counts = matrix( - ncol = nrow(object), dimnames = list(NULL, rownames(object)))), - reducedDims = list(RDA = object) - ) - return(object) -} - -# This function retrieves optional data that is used for creating an ellipses. -#' @importFrom scater retrieveCellInfo -.get_rda_ellipse_data <- function( - tse, reduced_dim, add.ellipse = TRUE, colour_by = color_by, - color_by = colour.by, colour.by = color.by, color.by = NULL, ...){ - # - if( !(add.ellipse %in% c(TRUE, FALSE, "fill", "color", "colour") && - length(add.ellipse) == 1L ) ){ - stop("'add.ellipse' must be one of c(TRUE, FALSE, 'fill', ", - "'color').", call. = FALSE) - } - if( !(is.null(colour_by) || .is_a_string(colour_by) && - colour_by %in% colnames(colData(tse)) ) ){ - stop("'colour_by' must be NULL or name of column from colData(x).", - call. = FALSE) - } - # - ellipse_data <- NULL - if( add.ellipse != FALSE && !is.null(colour_by) ){ - # Ellipse data is the same ordination data - ellipse_data <- as.data.frame(reduced_dim) - # Add sample metadata from colData - ellipse_data[[colour_by]] <- retrieveCellInfo(tse, colour_by)[["value"]] - attributes(ellipse_data)[["colour_by"]] <- colour_by - } - return(ellipse_data) -} - # This function retrieves data for creating vectors. Moreover, it wrangles the # vector data and controls what information is added to vector text or labels. .get_rda_vector_data <- function( @@ -507,43 +435,6 @@ setMethod("plotRDA", signature = c(x = "matrix"), return(vector_data) } -# This functions returns optional centroids for plotting. -.get_rda_centroids_data <- function( - reduced_dim, add.centroids = FALSE, ncomponents = 2L, ...){ - # - if( !.is_a_bool(add.centroids) ){ - stop("'add.centroids' must be TRUE or FALSE.", call. = FALSE) - } - if( !.is_an_integer(ncomponents) ){ - stop("'ncomponents' must be an integer.", call. = FALSE) - } - # - res <- if(add.centroids) .get_rda_attribute(reduced_dim, "centroids") - if( !is.null(res) ){ - res <- res[, seq_len(ncomponents), drop = FALSE] |> as.data.frame() - colnames(res) <- c("x", "y") - } - return(res) -} - -# This functions returns optional species scores for plotting. -.get_rda_species_data <- function( - reduced_dim, add.species = FALSE, ncomponents = 2L, ...){ - if( !.is_a_bool(add.species) ){ - stop("'add.species' must be TRUE or FALSE.", call. = FALSE) - } - if( !.is_an_integer(ncomponents) ){ - stop("'ncomponents' must be an integer.", call. = FALSE) - } - # - res <- if(add.species) .get_rda_attribute(reduced_dim, "species") - if( add.species ){ - res <- res[, seq_len(ncomponents), drop = FALSE] |> as.data.frame() - colnames(res) <- c("x", "y") - } - return(res) -} - # This function is used to fetch specified datatype from attributes if it # exists. .get_rda_attribute <- function(reduced_dim, attr_names){ @@ -556,137 +447,6 @@ setMethod("plotRDA", signature = c(x = "matrix"), return(res) } -# This function utilizes scater::plotReducedDim to create "baseplot". To where -# we can build the the plot. The idea is that the theme is similar in all -# ordination plots. -#' @importFrom scater plotReducedDim -#' @importFrom ggplot2 xlab ylab -.create_rda_baseplot <- function( - tse, dimred, reduced_dim, ncomponents = 2L, - add.expl.var = FALSE, expl.var = expl_var, expl_var = NULL, - colour_by = color_by, color_by = colour.by, - colour.by = color.by, color.by = NULL, ...){ - # - if( !.is_a_bool(add.expl.var) ){ - stop("'add.expl.var' must be TRUE or FALSE.", call. = FALSE) - } - if( !.is_an_integer(ncomponents) ){ - stop("'ncomponents' must be an integer.", call. = FALSE) - } - if( !( is.null(expl.var) || (is.numeric(expl.var) && - length(expl.var) == ncomponents )) ){ - stop("'expl.var' must be numeric vector with length ", ncomponents, - ".", call. = FALSE) - } - if( !(is.null(colour_by) || .is_a_string(colour_by) && - colour_by %in% colnames(colData(tse)) ) ){ - stop("'colour_by' must be NULL or name of column from colData(x).", - call. = FALSE) - } - # - # If specified, get explained variance - if( add.expl.var && is.null(expl.var) ){ - eigen_vals <- attr(reduced_dim, "eig") - # Convert to explained variance and take only first two components - expl_var <- eigen_vals / sum(eigen_vals) - expl_var <- expl_var[seq_len(ncomponents)]*100 - expl.var <- expl_var |> round(digits = 1) - } - # Create argument list - args <- c(list(object = tse, dimred = dimred, ncomponents = ncomponents, - colour_by = colour_by, percentVar = expl.var), list(...)) - # Remove additional arguments since plotReducedDim fails if we feed - # values that are not recognized - remove <- names(args) %in% c( - "add.significance", "add.expl.var", "add.ellipse", "add.vectors", - "vec.lab", "sep.group", "repl.underscore", "add.centroids", - "add.species", "ellipse.alpha", "ellipse.linewidth", "ellipse.linetype", - "confidence.level", "vec.size", "vec.color", "vec.colour", - "vec.linetype", "arrow.size", "min.segment.length", "label.color", - "label.colour", "label.size", "parse.labels", "vec.text", - "repel.labels", "position", "nudge_x", "nudge_y", "direction", - "max.overlaps", "check_overlap", "ignore.case") - args <- args[ !remove ] - # Get scatter plot with plotReducedDim --> keep theme similar between - # ordination methods - p <- do.call(plotReducedDim, args) - # Replace axis titles with names from result matrix - xlab <- colnames(reducedDim(tse, dimred))[[1L]] - ylab <- colnames(reducedDim(tse, dimred))[[2L]] - if( !is.null(expl.var) ){ - xlab <- paste0(xlab, " (", expl.var[[1L]], "%)") - ylab <- paste0(ylab, " (", expl.var[[2L]], "%)") - } - p <- p + xlab(xlab) + ylab(ylab) - return(p) -} - -# This function is used to create the plot. -.rda_plotter <- function(plot_data, ...){ - # Get the scatter plot - plot <- plot_data[["plot"]] - # Add ellipse - plot <- .rda_plotter_ellipse(plot, plot_data, ...) - # Add vectors - plot <- .rda_plotter_vector(plot, plot_data, ...) - # Add centroids - plot <- .rda_plotter_centroids_or_species(plot, plot_data, "centroids") - # Add species - plot <- .rda_plotter_centroids_or_species(plot, plot_data, "species_scores") - return(plot) -} - -# This function adds ellipse visualization. -.rda_plotter_ellipse <- function( - plot, plot_data, add.ellipse = TRUE, ellipse.alpha = 0.2, - ellipse.linewidth = 0.1, ellipse.linetype = 1, confidence.level = 0.95, - ...){ - # To disable "no visible binding for global variable" message in cmdcheck - color <- NULL - # - if( !(add.ellipse %in% c(TRUE, FALSE, "fill", "color", "colour") && - length(add.ellipse) == 1L ) ){ - stop("'add.ellipse' must be one of c(TRUE, FALSE, 'fill', ", - "'color').", call. = FALSE) - } - if( !.are_whole_numbers(ellipse.linetype) ){ - stop("'vec.linetype' must be a whole number.", call. = FALSE) - } - if ( !(is.numeric(ellipse.alpha) && ellipse.alpha > 0 && - ellipse.alpha < 1 ) ) { - stop("'ellipse.alpha' must be a number between 0 and 1.", call. = FALSE) - } - if ( !(is.numeric(ellipse.linewidth) && ellipse.linewidth > 0) ) { - stop("'ellipse.linewidth' must be a positive number.", call. = FALSE) - } - if( !(is.numeric(confidence.level) && confidence.level > 0 && - confidence.level < 1) ) { - stop("'confidence.level' must be a number between 0 and 1.", - call. = FALSE) - } - # - data <- plot_data[["ellipse_data"]] - if( !is.null(data) && length(data) > 2L){ - xvar <- colnames(data)[[1]] - yvar <- colnames(data)[[2]] - colour_var <- attributes(data)[["colour_by"]] - # Add ellipses to plot (fill or colour the edge) - fill <- add.ellipse %in% c(TRUE, "fill") - plot <- plot + stat_ellipse( - data = data, - mapping = aes( - x = .data[[xvar]], y = .data[[yvar]], - color = .data[[colour_var]], fill = after_scale(color)), - geom = "polygon", - linewidth = ellipse.linewidth, - linetype = ellipse.linetype, - level = confidence.level, - alpha = if(fill) ellipse.alpha else 0 - ) - } - return(plot) -} - # This function adds vector and text layer to the plot. #' @importFrom ggrepel geom_text_repel geom_label_repel .rda_plotter_vector <- function( @@ -696,7 +456,7 @@ setMethod("plotRDA", signature = c(x = "matrix"), vec.text = TRUE, repel.labels = TRUE, parse.labels = TRUE, add.significance = TRUE, vec.linetype = 1, min.segment.length = 5, position = NULL, nudge_x = NULL, nudge_y = NULL, direction = "both", - max.overlaps = 10, check_overlap = FALSE, ... + max.overlaps = 10, ... ){ # if ( !(is.numeric(vec.size) && vec.size > 0) ) { @@ -765,8 +525,6 @@ setMethod("plotRDA", signature = c(x = "matrix"), max.time = 0.5, max.iter = 10000, max.overlaps = max.overlaps, direction = direction, seed = NA, verbose = FALSE ) - } else if( !repel.labels && vec.text ){ - label_args <- c(label_args, check_overlap = check_overlap) } # Choose right function and call it FUN <- if( repel.labels && vec.text ) geom_text_repel @@ -777,19 +535,3 @@ setMethod("plotRDA", signature = c(x = "matrix"), } return(plot) } - -# This function adds centroids or species layer to the plot. -.rda_plotter_centroids_or_species <- function(plot, plot_data, type){ - # To disable "no visible binding for global variable" message in cmdcheck - x <- y <- NULL - data <- plot_data[[type]] - if( !is.null(data) ){ - plot <- plot + geom_point( - data, - mapping = aes(x = x, y = y), - shape = if(type == "centroids") 10L else 4L, - color = if(type == "centroids") "blue" else "red", - ) - } - return(plot) -} diff --git a/R/plotJointRPCA.R b/R/plotJointRPCA.R new file mode 100644 index 00000000..3611bd19 --- /dev/null +++ b/R/plotJointRPCA.R @@ -0,0 +1,339 @@ +#' @name +#' plotJointRPCA +#' +#' @title +#' Visualize Joint-RPCA results +#' +#' @description +#' Creates a two-dimensional ordination plot from Joint-RPCA results stored in a +#' \code{SingleCellExperiment} or \code{MultiAssayExperiment}. In addition to the +#' sample ordination, feature loadings can be visualized as vectors. +#' +#' @details +#' This function is a wrapper around \code{\link[=plotOrdination]{plotOrdination}} +#' for Joint-RPCA results. Consequently, most graphical parameters are passed +#' directly to \code{plotOrdination()}, including options for colouring, +#' grouping, faceting and adding ellipses. See +#' \code{\link[=plotOrdination]{plotOrdination}} for a complete description of +#' these arguments. +#' +#' Feature loadings are plotted as vectors originating from the origin. By +#' default, only the longest loading vectors are shown for each data layer to +#' improve readability. +#' +#' @return +#' A \code{ggplot2} object. +#' +#' @param x A +#' \code{\link[SummarizedExperiment:SummarizedExperiment-class]{SummarizedExperiment}} +#' or +#' \code{\link[MultiAssayExperiment:MultiAssayExperiment-class]{MultiAssayExperiment}} +#' object containing Joint-RPCA results. +#' +#' @param dimred \code{character(1)}: Specifies the name of the Joint-RPCA +#' result to plot. +#' +#' @param add.vectors \code{logical(1)} or \code{character}. If +#' \code{TRUE}, feature loading vectors are added. Alternatively, a character +#' vector can be supplied to display only features whose names match the given +#' pattern(s). (Default: \code{TRUE}) +#' +#' @param ... +#' Additional arguments passed to +#' \code{\link[=plotOrdination]{plotOrdination}} and to the feature-vector +#' plotting. +#' Commonly used Joint-RPCA-specific arguments include: +#' \itemize{ +#' \item \code{ntop}: \code{integer(1)} or \code{NULL}. Maximum number of +#' loading vectors shown per data layer. The vectors with the largest lengths +#' are retained. Set to \code{NULL} to display all vectors. (Default: \code{10}) +#' } +#' +#' @examples +#' data("HintikkaXOData") +#' mae <- HintikkaXOData +#' +#' mae[[1]] <- transformAssay( +#' mae[[1]], +#' assay.type = "counts", +#' method = "rclr", +#' impute = FALSE +#' ) +#' mae[[2]] <- transformAssay( +#' mae[[2]], +#' assay.type = "nmr", +#' method = "log10" +#' ) +#' +#' mae <- addJointRPCA( +#' mae, +#' experiments = c(1, 2), +#' assay.types = c("rclr", "log10") +#' ) +#' +#' # Basic Joint-RPCA plot +#' plotJointRPCA(mae, "JointRPCA") +#' +#' # Colour samples by metadata +#' plotJointRPCA(mae, "JointRPCA", colour.by = "Fat") +#' +#' # Add confidence ellipses +#' plotJointRPCA( +#' mae, +#' "JointRPCA", +#' colour.by = "Fat", +#' add.ellipse = TRUE +#' ) +#' +#' # Show all loading vectors +#' plotJointRPCA( +#' mae, +#' "JointRPCA", +#' ntop = NULL +#' ) +#' +#' # Display only selected feature vectors +#' plotJointRPCA( +#' mae, +#' "JointRPCA", +#' add.vectors = c("Bacteroides", "Roseburia") +#' ) +#' +#' # Use boxed labels without repelling +#' plotJointRPCA( +#' mae, +#' "JointRPCA", +#' text.labels = FALSE, +#' repel.labels = FALSE +#' ) +#' +#' @seealso +#' \itemize{ +#' \item \code{\link[=plotOrdination]{plotOrdination}} +#' } +#' +NULL + +#' @rdname plotJointRPCA +#' @export +setMethod("plotJointRPCA", signature = c(x = "MultiAssayExperiment"), + function(x, dimred, ...){ + .check_metadata_present(dimred, x) + # Construct TreeSE from results so that we can use TreeSE function + # (and thus plotOrdination) + res <- metadata(x)[[dimred]] + # + col_data <- x |> colData() + col_data <- col_data[ + match(rownames(res), rownames(col_data)), , drop = FALSE] + rownames(col_data) <- rownames(res) + # + tse <- TreeSummarizedExperiment( + assays = SimpleList(counts = matrix( + ncol = nrow(res), dimnames = list(NULL, rownames(res)))), + colData = col_data, + reducedDims = list(JointRPCA = res) + ) + # Plot + p <- plotJointRPCA(tse, dimred = "JointRPCA", ...) + return(p) + } +) + +#' @rdname plotJointRPCA +#' @export +setMethod("plotJointRPCA", signature = c(x = "SingleCellExperiment"), + function(x, dimred, add.vectors = TRUE, ...){ + .check_dimred_present(dimred, x) + if( !inherits(reducedDim(x, dimred), "JointRPCA") ){ + stop("The plotted results must be a class 'JointRPCA'.", + call. = FALSE) + } + if( !( .is_a_bool(add.vectors) || is.character(add.vectors)) ){ + stop("'add.vectors must be TRUE or FALSE or character vector.", + call. = FALSE) + } + p <- plotOrdination(x, dimred, ...) + vector_data <- .get_joint_rpca_vector_data(x, dimred, add.vectors, ...) + if( !is.null(vector_data ) ){ + p <- .add_ordination_vectors(p, vector_data, ...) + } + return(p) + } +) + +################################ HELP FUNCTIONS ################################ + +# Get feature loadings that will be plotted as vectors +#' @importFrom dplyr group_by slice_max ungroup +.get_joint_rpca_vector_data <- function( + tse, reduced_dim, add.vectors, ignore.case = FALSE, ntop = 10, ...){ + if( !.is_a_bool(ignore.case) ){ + stop("'ignore.case' must be TRUE or FALSE.", call. = FALSE) + } + if( !(.is_an_integer(ntop) || is.null(ntop)) ){ + stop("'ntop' must be an integer.", call. = FALSE) + } + # + # Get vector data, i.e, rotations or loadings + res <- reducedDim(tse, reduced_dim) + vector_data <- if(!(.is_a_bool(add.vectors) && !add.vectors)) + attributes(res)[["rotation"]] + # If user wanted to plot them, wrangle the data + if( !is.null(vector_data) ){ + # Add labels and layer name + vector_data <- as.data.frame(vector_data) + vector_data[["vector_label"]] <- rownames(vector_data) + n_features <- attr(res, "n_features") + vector_data[["Layer"]] <- rep( + names(n_features), + times = unname(n_features) + ) + + # Subset vectors by selecting only those ones that user has specified + if( is.character(add.vectors) ){ + add.vectors <- paste0(add.vectors, collapse = "|") + keep <- vapply(rownames(vector_data), function(x) + grepl(add.vectors, x, perl = TRUE, ignore.case = ignore.case), + logical(1L)) + vector_data <- vector_data[keep, ] + } + + # Keep only the n longest vectors per layer. The idea is to show only + # features that are highly associated with the showed ordination space. + if (!is.null(ntop) && nrow(vector_data) > 0L) { + xvar <- colnames(vector_data)[1] + yvar <- colnames(vector_data)[2] + + vector_data[["length_temporary"]] <- + sqrt(vector_data[[xvar]]^2 + vector_data[[yvar]]^2) + + vector_data <- vector_data |> + group_by(Layer) |> + slice_max(length_temporary, n = ntop, with_ties = FALSE) |> + ungroup() + + vector_data[["length_temporary"]] <- NULL + } + + # If all vectors were removed, give NULL + if( nrow(vector_data) == 0L ){ + vector_data <- NULL + } + } + + return(vector_data) +} + +# Add feature loadings as vectors +#' @importFrom ggrepel geom_text_repel geom_label_repel +#' @importFrom grid arrow unit +.add_ordination_vectors <- function( + p, + vector_data, + vec.size = 0.5, + arrow.size = 0.25, + label.size = 4, + vec.color = "black", + label.color = "black", + vec.text = TRUE, + text.labels = TRUE, + repel.labels = TRUE, + min.segment.length = 0.5, + box.padding = 0.25, + point.padding = 1e-06, + force = 1, + force_pull = 1, + max.time = 0.5, + max.iter = 10000, + max.overlaps = 10, + direction = "both", + seed = NA, + ... +){ + if ( !(is.numeric(vec.size) && vec.size > 0) ) { + stop("'vec.size' must be a positive number.", call. = FALSE) + } + if ( !(is.numeric(arrow.size) && arrow.size > 0) ) { + stop("'arrow.size' must be a positive number.", call. = FALSE) + } + if ( !(is.numeric(label.size) && label.size > 0) ) { + stop("'label.size' must be a positive number.", call. = FALSE) + } + if ( !.is_non_empty_string(vec.color) ) { + stop("'vec.color' must be a non-empty string specifying a colour", + call. = FALSE) + } + if ( !.is_non_empty_string(label.color) ) { + stop("'label.color' must be a non-empty string specifying a colour", + call. = FALSE) + } + if( !.is_a_bool(vec.text) ){ + stop("'vec.text' must be TRUE or FALSE.", call. = FALSE) + } + if( !.is_a_bool(repel.labels) ){ + stop("'repel.labels' must be TRUE or FALSE.", call. = FALSE) + } + # + # Names of the ordination axes + xvar <- colnames(vector_data)[1] + yvar <- colnames(vector_data)[2] + + # Draw vectors from the origin to the feature coordinates + p <- p + + geom_segment( + data = vector_data, + aes( + x = 0, + y = 0, + xend = .data[[xvar]], + yend = .data[[yvar]], + linetype = .data[["Layer"]] + ), + arrow = arrow(length = unit(arrow.size, "cm")), + linewidth = vec.size, + colour = vec.color, + inherit.aes = FALSE + ) + + # Select the appropriate text/label geometry + FUN <- if (repel.labels) { + if (text.labels) geom_text_repel else geom_label_repel + } else { + if (text.labels) geom_text else geom_label + } + + # Common arguments shared by all text geometries + label_args <- list( + data = vector_data, + mapping = aes( + x = .data[[xvar]], + y = .data[[yvar]], + label = .data[["vector_label"]] + ), + colour = label.color, + size = label.size, + inherit.aes = FALSE + ) + + # Add geometry-specific arguments + if (repel.labels) { + label_args <- c(label_args, list( + min.segment.length = min.segment.length, + box.padding = box.padding, + point.padding = point.padding, + force = force, + force_pull = force_pull, + max.time = max.time, + max.iter = max.iter, + max.overlaps = max.overlaps, + direction = direction, + seed = seed + )) + } + + # Add feature labels + p <- p + do.call(FUN, label_args) + + return(p) +} diff --git a/R/plotOrdination.R b/R/plotOrdination.R new file mode 100644 index 00000000..11ad2d0b --- /dev/null +++ b/R/plotOrdination.R @@ -0,0 +1,831 @@ +#' @name +#' plotOrdination +#' +#' @title +#' Visualize ordination results +#' +#' @description +#' Creates a two-dimensional ordination plot from a reduced dimension result +#' stored in a +#' \code{\link[SingleCellExperiment:reducedDims]{SingleCellExperiment}} +#' object. Samples can be coloured, filled, shaped, sized, grouped or faceted +#' using sample metadata or feature abundances, and additional graphical +#' elements such as confidence ellipses, centroids, vectors and density +#' estimates can be added. +#' +#' @details +#' This function provides a unified interface for visualizing ordination methods +#' such as PCA, PCoA, MDS, t-SNE, UMAP, RDA and CCA. The plotted coordinates are +#' retrieved from a reduced dimension result stored in +#' \code{reducedDim(x)}. +#' +#' +#' @return +#' A \code{ggplot2} object. +#' +#' @param x a +#' \code{\link[SummarizedExperiment:SummarizedExperiment-class]{SummarizedExperiment}} +#' object. +#' +#' @param dimred \code{character scalar}. Name of the reduced dimension result +#' stored in \code{reducedDim(x)} to visualize. +#' +#' @param ... Additional parameters for plotting. +#' \itemize{ +#' \item \code{ncomponents}: \code{integer vector} of length 2 or +#' \code{integer scalar}. Specifies which ordination components are plotted. +#' If a scalar is provided, the first two components are used. +#' (Default: \code{2L}) +#' +#' \item \code{colour.by}: \code{NULL} or \code{character scalar}. Specifies a +#' variable from \code{colData(x)} or a feature from \code{rownames(x)} used +#' to colour observations. Feature abundances are taken from +#' \code{assay.type}. (Default: \code{NULL}) +#' +#' \item \code{fill.by}: \code{NULL} or \code{character scalar}. Specifies a +#' variable from \code{colData(x)} or a feature from \code{rownames(x)} used +#' to fill observations or ellipses. Feature abundances are taken from +#' \code{assay.type}. Cannot be used together with +#' \code{add.density = TRUE}. (Default: \code{NULL}) +#' +#' \item \code{shape.by}: \code{NULL} or \code{character scalar}. Specifies a +#' categorical variable from \code{colData(x)} used for point shapes. +#' (Default: \code{NULL}) +#' +#' \item \code{size.by}: \code{NULL} or \code{character scalar}. Specifies a +#' variable from \code{colData(x)} used for point sizes. +#' (Default: \code{NULL}) +#' +#' \item \code{group.by}: \code{NULL} or \code{character scalar}. Specifies a +#' categorical variable from \code{colData(x)} used for grouping when drawing +#' ellipses, centroids and centroid vectors. (Default: \code{NULL}) +#' +#' \item \code{linetype.by}: \code{NULL} or \code{character scalar}. +#' Specifies a categorical variable from \code{colData(x)} used for ellipse +#' line types. (Default: \code{NULL}) +#' +#' \item \code{pair.by}: \code{NULL} or \code{character scalar}. Specifies a +#' variable from \code{colData(x)} identifying observations that should be +#' connected by lines. (Default: \code{NULL}) +#' +#' \item \code{sort.by}: \code{NULL} or \code{character scalar}. Specifies a +#' variable from \code{colData(x)} used to order observations before drawing +#' connecting lines. (Default: \code{NULL}) +#' +#' \item \code{facet.by}: \code{NULL} or \code{character scalar}. Specifies a +#' categorical variable from \code{colData(x)} used to split the plot into +#' facets. (Default: \code{NULL}) +#' +#' \item \code{assay.type}: \code{character scalar}. Name of the assay used +#' when \code{colour.by} or \code{fill.by} specifies a feature. +#' (Default: \code{"counts"}) +#' +#' \item \code{add.points}: \code{logical scalar}. Whether to draw sample +#' points. (Default: \code{TRUE}) +#' +#' \item \code{add.ellipse}: \code{logical scalar}. Whether to draw confidence +#' ellipses around groups. (Default: \code{FALSE}) +#' +#' \item \code{add.density}: \code{logical scalar}. Whether to draw a +#' two-dimensional density estimate in the background. +#' (Default: \code{FALSE}) +#' +#' \item \code{add.centroids}: \code{logical scalar}. Whether to draw group +#' centroids. (Default: \code{FALSE}) +#' +#' \item \code{add.centroids.lines}: \code{logical scalar}. Whether to connect +#' observations to their group centroids. (Default: \code{FALSE}) +#' +#' \item \code{add.vectors}: \code{logical scalar}. Whether to draw vectors +#' from the global centroid to group centroids. +#' (Default: \code{FALSE}) +#' +#' \item \code{add.rotation}: \code{logical scalar}. Whether to draw rotation +#' (species score) coordinates if available in the ordination result. +#' (Default: \code{add.species}) +#' +#' \item \code{add.species}: \code{logical scalar}. Alias for +#' \code{add.rotation}. (Default: \code{FALSE}) +#' +#' \item \code{add.expl.var}: \code{logical scalar}. Whether to append the +#' percentage of explained variance to the axis labels when available. +#' (Default: \code{FALSE}) +#' +#' \item \code{scales}: \code{character scalar}. Scaling used for faceted +#' plots. Passed to \code{ggplot2::facet_wrap()}. (Default: +#' \code{"fixed"}) +#' +#' \item \code{xlab}, \code{ylab}: \code{character scalar}. Axis labels. +#' Defaults to the ordination component names. +#' +#' \item \code{panel.by.eigen}: \code{logical scalar}. Whether to scale the +#' panel aspect ratio according to the eigenvalues of the ordination when +#' available. (Default: \code{TRUE}) +#' +#' \item \code{point.shape}: Shape used for points. +#' (Default: \code{19}) +#' +#' \item \code{point.alpha}: \code{numeric scalar}. Transparency of points. +#' Must be between 0 and 1. (Default: \code{0.4}) +#' +#' \item \code{ellipse.alpha}: \code{numeric scalar}. Transparency of ellipse +#' fills. Must be between 0 and 1. (Default: \code{0.2}) +#' +#' \item \code{ellipse.linewidth}: \code{numeric scalar}. Line width of +#' ellipse borders. (Default: \code{0.5} or \code{0} when +#' \code{fill.by} is specified.) +#' +#' \item \code{ellipse.linetype}: Integer specifying the ellipse line type. +#' (Default: \code{1}) +#' +#' \item \code{confidence.level}: \code{numeric scalar}. Confidence level used +#' for ellipse calculation. Must be between 0 and 1. +#' (Default: \code{0.95}) +#' +#' \item \code{adjust}: \code{numeric scalar}. Multiplicative adjustment for +#' the bandwidth used in the background density estimate. +#' (Default: \code{1}) +#' } +#' +#' @examples +#' data("Tito2024QMP") +#' tse <- Tito2024QMP +#' +#' # Compute relative abundances and an MDS ordination +#' tse <- transformAssay(tse, method = "relabundance") +#' tse <- addMDS(tse, assay.type = "relabundance", method = "bray", ncomponents = 50) +#' +#' # Basic ordination plot +#' plotOrdination(tse, "MDS") +#' +#' # Colour samples by a sample-level variable +#' plotOrdination(tse, "MDS", colour.by = "diagnosis") +#' +#' # Colour samples by the abundance of a single feature +#' plotOrdination( +#' tse, "MDS", +#' colour.by = rownames(tse)[1], +#' assay.type = "relabundance") +#' +#' # Use multiple aesthetics simultaneously +#' plotOrdination( +#' tse, "MDS", +#' colour.by = "diagnosis", +#' shape.by = "colonoscopy" +#' ) +#' +#' # Add confidence ellipses +#' plotOrdination( +#' tse, "MDS", +#' colour.by = "diagnosis", +#' group.by = "diagnosis", +#' add.ellipse = TRUE +#' ) +#' +#' # Fill ellipses instead of colouring them +#' plotOrdination( +#' tse, "MDS", +#' fill.by = "diagnosis", +#' group.by = "diagnosis", +#' add.ellipse = TRUE +#' ) +#' +#' # Show explained variance in the axis labels +#' plotOrdination( +#' tse, "MDS", +#' colour.by = "diagnosis", +#' add.expl.var = TRUE +#' ) +#' +#' # Add group centroids +#' plotOrdination( +#' tse, "MDS", +#' colour.by = "diagnosis", +#' group.by = "diagnosis", +#' add.centroids = TRUE +#' ) +#' +#' # Connect samples to their group centroids +#' plotOrdination( +#' tse, "MDS", +#' colour.by = "diagnosis", +#' group.by = "diagnosis", +#' add.centroids.lines = TRUE +#' ) +#' +#' # Show vectors from the global centroid to each group centroid +#' plotOrdination( +#' tse, "MDS", +#' colour.by = "diagnosis", +#' group.by = "diagnosis", +#' add.vectors = TRUE +#' ) +#' +#' # Draw a background density estimate +#' plotOrdination( +#' tse, "MDS", +#' add.density = TRUE +#' ) +#' +#' # Split the plot into facets +#' plotOrdination( +#' tse, "MDS", +#' colour.by = "diagnosis", +#' facet.by = "colonoscopy" +#' ) +#' +#' # Plot different ordination components +#' plotOrdination( +#' tse, "MDS", +#' ncomponents = c(2, 3) +#' ) +#' +#' # Customize point appearance +#' plotOrdination( +#' tse, "MDS", +#' colour.by = "diagnosis", +#' point.shape = 17, +#' point.alpha = 0.8 +#' ) +#' +#' @seealso +#' \itemize{ +#' \item \code{\link[scater:plotReducedDim]{scater::plotReducedDim}} +#' \item \code{\link[=plotCCA]{plotCCA}} +#' } +#' +NULL + +#' @rdname plotOrdination +#' @export +setMethod("plotOrdination", signature = c(x = "SingleCellExperiment"), + function(x, dimred, ...){ + args <- .check_ordination_input(x, dimred, ...) + df <- do.call(.get_ordination_data, args) + p <- .ordination_plotter(df, ...) + return(p) + } +) + +################################ HELP FUNCTIONS ################################ + +# This method checks that the input is in correct format +.check_ordination_input <- function( + x, dimred, + ncomponents = 2L, + colour.by = color.by, color.by = NULL, + fill.by = NULL, + shape.by = NULL, + size.by = NULL, + group.by = NULL, + linetype.by = NULL, + pair.by = NULL, + sort.by = NULL, + facet.by = NULL, + assay.type = "counts", + add.points = TRUE, add.ellipse = FALSE, add.density = FALSE, + add.centroids = FALSE, add.centroids.lines = FALSE, add.vectors = FALSE, + add.rotation = add.species, add.species = FALSE, add.expl.var = FALSE, + ...){ + # Check that dimred can be found + .check_dimred_present(dimred, x) + # Check that ncomponents is correct. We can only visualize 2 components. + if( .is_an_integer(ncomponents) ){ + ncomponents <- seq_len(ncomponents) + } + if( !(.is_integer(ncomponents) && length(ncomponents) == 2L && + all(ncomponents > 0L & ncomponents <= ncol(reducedDim(x, dimred)))) + ){ + stop("'ncomponents' must specify columns from reducedDim(x, dimred) ", + "with integer values.", call. = FALSE) + } + + # Check aesthetic variables + temp <- .check_metadata_variable(x, colour.by, FALSE, TRUE, FALSE, TRUE) + temp <- .check_metadata_variable(x, fill.by, FALSE, TRUE, FALSE, TRUE) + temp <- .check_metadata_variable(x, shape.by, FALSE, TRUE, FALSE, FALSE) + temp <- .check_metadata_variable(x, size.by, FALSE, TRUE, FALSE, FALSE) + temp <- .check_metadata_variable(x, linetype.by, FALSE, TRUE, FALSE, FALSE) + temp <- .check_metadata_variable(x, group.by, FALSE, TRUE, FALSE, FALSE) + temp <- .check_metadata_variable(x, pair.by, FALSE, TRUE, FALSE, FALSE) + temp <- .check_metadata_variable(x, sort.by, FALSE, TRUE, FALSE, FALSE) + temp <- .check_metadata_variable(x, facet.by, FALSE, TRUE, FALSE, FALSE) + # If colour.by specifies rowname, we check assay.type as the abundance + # values are used for coloring + if( !is.null(colour.by) && colour.by %in% rownames(x) ){ + temp <- .check_assay_present(assay.type, x) + } + # Check other flags + if( !.is_a_bool(add.points) ){ + stop("'add.points' must be TRUE or FALSE.", call. = FALSE) + } + if( !.is_a_bool(add.ellipse) ){ + stop("'add.ellipse' must be TRUE or FALSE.", call. = FALSE) + } + if( !.is_a_bool(add.density) ){ + stop("'add.density' must be TRUE or FALSE.", call. = FALSE) + } + if( !.is_a_bool(add.centroids) ){ + stop("'add.centroids' must be TRUE or FALSE.", call. = FALSE) + } + if( !.is_a_bool(add.centroids.lines) ){ + stop("'add.centroids.lines' must be TRUE or FALSE.", call. = FALSE) + } + if( !.is_a_bool(add.vectors) ){ + stop("'add.vectors' must be TRUE or FALSE.", call. = FALSE) + } + if( !.is_a_bool(add.rotation) ){ + stop("'add.rotation' must be TRUE or FALSE.", call. = FALSE) + } + if( !.is_a_bool(add.expl.var) ){ + stop("'add.expl.var' must be TRUE or FALSE.", call. = FALSE) + } + + # add.density cannot be specified simultaneously with fill.by as they both + # are using fill aesthetic and we can have only on fill scale. + if( add.density && !is.null(fill.by) ){ + stop("Both 'add.density' and 'fill.by' cannot be specified ", + "simultaneously.", call. = FALSE) + } + # Put all the arguments into list that be fed to the data retrieval function + args <- list( + x = x, dimred = dimred, + ncomponents = ncomponents, + colour.by = colour.by, + fill.by = fill.by, + shape.by = shape.by, + size.by = size.by, + group.by = group.by, + linetype.by = linetype.by, + pair.by = pair.by, + sort.by = sort.by, + facet.by = facet.by, + assay.type = assay.type, + add.expl.var = add.expl.var + ) + args <- c(args, list(...)) + return(args) +} + +# This function retrieves the data from reducedDim +.get_ordination_data <- function( + x, dimred, ncomponents, colour.by, fill.by, shape.by, size.by, group.by, + linetype.by, pair.by, sort.by, facet.by, assay.type, + add.expl.var = FALSE, ...){ + # Get data and store the original attributes that might include rotation + # data, for instance + df <- reducedDim(x, dimred) + orig_attributes <- attributes(df) + orig_attributes <- orig_attributes[ + !names(orig_attributes) %in% c("dim", "dimnames") ] + # Subset the data to include specified columns to plot + df <- df[, ncomponents] + # Add colnames if they are not present + if( is.null(colnames(df)) ){ + colnames(df) <- paste0(dimred, ncomponents) + } + df <- df |> as.data.frame() + # Take names of the columns that will be plotted + x_var <- colnames(df)[[1L]] + y_var <- colnames(df)[[2L]] + + # Get rotation data and put it in correct format + rotation <- NULL + rotation_names <- c("rotation", "species") + if( any(rotation_names %in% names(orig_attributes)) ){ + rotation_names <- rotation_names[ + rotation_names %in% names(orig_attributes)][[1L]] + rotation <- orig_attributes[[rotation_names]] |> as.data.frame() + rotation <- rotation[, ncomponents, drop = FALSE] + colnames(rotation) <- colnames(df) + } + + xlab <- x_var + ylab <- y_var + eigen_names <- c("eig", "percentVar") + if( add.expl.var && any(eigen_names %in% names(orig_attributes)) ){ + eigen_names <- eigen_names[ + eigen_names %in% names(orig_attributes)][[1L]] + eigen <- orig_attributes[[eigen_names]] + if( eigen_names %in% c("eig") ){ + eigen <- eigen * 100 + } + xlab <- paste0( + xlab, " (", + round(eigen[ncomponents][[1L]], 1), + "%)" + ) + ylab <- paste0( + ylab, " (", + round(eigen[ncomponents][[2L]], 1), + "%)" + ) + } else if( add.expl.var ){ + warning("No explained variance found from the data.", call. = FALSE) + } + + # List data that is fetched from colData + cols <- c( + shape.by, size.by, group.by, linetype.by, pair.by, sort.by, facet.by) + # colour.by and fill.by can also specify a feature and its abundance. + # Either get their column name to fetch from colData or directly add them + # to data. + if( !is.null(colour.by) && colour.by %in% colnames(colData(x)) ){ + cols <- c(cols, colour.by) + } else if( !is.null(colour.by) ){ + df[[colour.by]] <- assay(x, assay.type)[colour.by, ] + } + if( !is.null(fill.by) && fill.by %in% colnames(colData(x)) ){ + cols <- c(cols, fill.by) + } else if( !is.null(fill.by) ){ + df[[fill.by]] <- assay(x, assay.type)[fill.by, ] + } + # Get colData and merge it with the data + cols <- cols |> unique() + cd <- colData(x)[, cols, drop = FALSE] + df <- cbind(df, cd) + + # Check that the values are correct + if( !is.null(group.by) && is.numeric(df[[group.by]]) ){ + stop("Values specified by 'group.by' must be categorical.", + call. = FALSE) + } + if( !is.null(linetype.by) && is.numeric(df[[linetype.by]]) ){ + stop("Values specified by 'linetype.by' must be categorical.", + call. = FALSE) + } + if( !is.null(facet.by) && is.numeric(df[[facet.by]]) ){ + stop("Values specified by 'facet.by' must be categorical.", + call. = FALSE) + } + + # Sort the data. For instance, if we want to add arrows between consecutive + # time points, this is essential step. + if( !is.null(sort.by) ){ + df <- df[order(df[[sort.by]]), , drop = FALSE] + } + + # Calculate centroid data + grouping_var <- c( + colour.by, facet.by, group.by, fill.by) |> unique() + df_centroids <- df |> + group_by(across(all_of(grouping_var))) |> + summarise( + x_centroid = mean(.data[[x_var]], na.rm = TRUE), + y_centroid = mean(.data[[y_var]], na.rm = TRUE), + .groups = "drop" + ) + # Add global mean + df_centroids[["x_global"]] <- mean(df[[x_var]], na.rm = TRUE) + df_centroids[["y_global"]] <- mean(df[[y_var]], na.rm = TRUE) + + # Add additional information to attributes of df. + attributes(df) <- c( + attributes(df), + x = x_var, + y = y_var, + colour.by = colour.by, + fill.by = fill.by, + shape.by = shape.by, + size.by = size.by, + group.by = group.by, + linetype.by = linetype.by, + pair.by = pair.by, + sort.by = sort.by, + facet.by = facet.by, + assay.type = assay.type, + xlab = xlab, + ylab = ylab + ) + attr(df, "rotation") <- rotation + attr(df, "centroids") <- df_centroids + attr(df, "eigen") <- eigen + return(df) +} + +# This method is the main plotter function. +.ordination_plotter <- function( + df, scales = "fixed", add.points = TRUE, add.ellipse = FALSE, + add.density = FALSE, add.centroids = FALSE, add.centroids.lines = FALSE, + add.vectors = FALSE, add.rotation = add.species, add.species = FALSE, + ...){ + # Initialize the plot + p <- ggplot(df, aes( + x = .data[[attributes(df)[["x"]]]], + y = .data[[attributes(df)[["y"]]]], + colour = if( !is.null(attributes(df)[["colour.by"]]) ) + .data[[attributes(df)[["colour.by"]]]] + )) + grouping_var <- c( + attributes(df)[["colour.by"]], attributes(df)[["facet.by"]], + attributes(df)[["group.by"]], attributes(df)[["fill.by"]]) |> unique() + # Add points, i.e., samples + if( add.points ){ + p <- .add_ordination_points(p, df, ...) + } + # Add ellipses + if( add.ellipse && length(grouping_var) > 0L ){ + p <- .add_ellipse_to_ordination(p, df, ...) + } + # Connect samples of subjects, groups etc + if( !is.null(attributes(df)[["pair.by"]]) ){ + p <- .connect_ordination_points(p, df) + } + # Add density to background + if( add.density ){ + # Background density + p <- .add_point_density(p, df) + } + # Add group centroids + if( add.centroids && length(grouping_var) > 0L ){ + p <- .add_centroids(p, df, ...) + } + # Add lines connecting points and group centroids + if( add.centroids.lines && length(grouping_var) > 0L ){ + p <- .add_centroid_lines(p, df, grouping_var, ...) + } + # Add species scores + if( add.rotation ){ + p <- .add_rotation(p, df) + } + # Add vectors from global centroid to group centroids + if( add.vectors && length(grouping_var) > 0L ){ + p <- .add_centroids_vector(p, df, grouping_var, ...) + } + # If facetting was specified, split plot to separate panels + if( !is.null(attributes(df)[["facet.by"]]) ){ + p <- p + + facet_wrap( + ~ .data[[attributes(df)[["facet.by"]]]], + scales = scales + ) + } + # Adjust theme + p <- .adjust_ordination_theme(p, df, ...) + return(p) +} + +# Add points for ordination plot +.add_ordination_points <- function( + p, df, point.shape = 19, point.alpha = 0.4, ...){ + p <- .add_points_layer( + p, df, + x = attributes(df)[["x"]], + y = attributes(df)[["y"]], + point.shape = point.shape, + point.alpha = point.alpha, + ...) + return(p) +} + +# This function adds ellipse visualization. +.add_ellipse_to_ordination <- function( + p, df, + ellipse.alpha = 0.2, + ellipse.linewidth = if(is.null(attributes(df)[["fill.by"]])) 0.5 else 0, + ellipse.linetype = 1, + confidence.level = 0.95, + ...){ + # To disable "no visible binding for global variable" message in cmdcheck + color <- NULL + # + if( !.are_whole_numbers(ellipse.linetype) ){ + stop("'vec.linetype' must be a whole number.", call. = FALSE) + } + if ( !(is.numeric(ellipse.alpha) && ellipse.alpha > 0 && + ellipse.alpha < 1 ) ) { + stop("'ellipse.alpha' must be a number between 0 and 1.", call. = FALSE) + } + if ( !(is.numeric(ellipse.linewidth) && ellipse.linewidth >= 0) ) { + stop("'ellipse.linewidth' must be a positive number.", call. = FALSE) + } + if( !(is.numeric(confidence.level) && confidence.level > 0 && + confidence.level < 1) ) { + stop("'confidence.level' must be a number between 0 and 1.", + call. = FALSE) + } + # + # Get all the arguments. User can fill and colour the ellipses separately. + # However, in most of the cases that might not make sense, but it is still + # made possible. + args <- list( + mapping = aes( + group = if( !is.null(attributes(df)[["group.by"]]) ) + .data[[attributes(df)[["group.by"]]]], + colour = if( !is.null(attributes(df)[["colour.by"]]) ) + .data[[attributes(df)[["colour.by"]]]], + fill = if( !is.null(attributes(df)[["fill.by"]]) ) + .data[[attributes(df)[["fill.by"]]]], + linetype = if( !is.null(attributes(df)[["linetype.by"]]) ) + .data[[attributes(df)[["linetype.by"]]]], + ), + geom = "polygon", + linewidth = ellipse.linewidth, + level = confidence.level, + alpha = if( is.null(attributes(df)[["fill.by"]]) ) 0 else ellipse.alpha + ) + # If user did not specify coloring, add black border to ellipses + if( is.null(attributes(df)[["colour.by"]]) ){ + + } + if( is.null(attributes(df)[["linetype.by"]]) ){ + args[["linetype"]] = ellipse.linetype + } + p <- p + do.call(stat_ellipse, args) + return(p) +} + +# This method connects points with a line or arrow. +.connect_ordination_points <- function(p, df){ + # Get arguments. If user sorted the data, use directed arrow. + args <- list( + mapping = aes( + x = .data[[attributes(df)[["x"]]]], + y = .data[[attributes(df)[["y"]]]], + group = .data[[attributes(df)[["pair.by"]]]], + ), + arrow = if(!is.null(attributes(df)[["sort.by"]])) + arrow(length = unit(0.2, "cm"), type = "closed"), + alpha = 0.4 + ) + args <- args[ lengths(args) > 0 ] + p <- p + do.call(geom_path, args) + return(p) +} + +# This method adds group +.add_centroids <- function(p, df, ...){ + df_centroids <- attributes(df)[["centroids"]] + # Add centroids + p <- p + + geom_point( + data = df_centroids, + aes(x = x_centroid, y = y_centroid), + size = 5, + shape = 4, # cross + stroke = 2 # thick border + ) + return(p) +} + +# This method adds group +.add_centroid_lines <- function(p, df, grouping_var, ...){ + df_centroids <- attributes(df)[["centroids"]] + # Add centroids data to original df + df <- df %>% + dplyr::left_join(df_centroids, by = grouping_var) + # Connect points with centroids + p <- p + geom_segment(data = df, aes( + x = x_centroid, y = y_centroid, + xend = .data[[attributes(df)[["x"]]]], + yend = .data[[attributes(df)[["y"]]]] + ), alpha = 0.4) + return(p) +} + +# This methods creates vectors that start from global mean and ends to group +# centroids. This shows how the covariate correlates with the ordination. +#' @importFrom ggrepel geom_label_repel +.add_centroids_vector <- function(p, df, grouping_var, ...){ + df_centroids <- attributes(df)[["centroids"]] + # Visualize vectors + p <- p + geom_segment( + data = df_centroids, + aes( + x = x_global, y = y_global, + xend = x_centroid, yend = y_centroid, + ), + arrow = arrow(length = unit(0.2, "cm"), type = "closed"), + size = 1 + ) + # Add label to denote which vector belongs to which group + p <- p + geom_label_repel( + data = df_centroids, + mapping = aes( + x = x_centroid, + y = y_centroid, + label = .data[[grouping_var]]) + ) + return(p) +} + +# This methods adds "species scores", i.e., the coordinates of species in the +# ordination. +.add_rotation <- function(p, df){ + # The points are added only if the rotation is present in the data + if( !is.null(attributes(df)[["rotation"]]) ){ + df_species <- attributes(df)[["rotation"]] + # Add points + p <- p + geom_point( + data = df_species, + mapping = aes( + x = .data[[attributes(df)[["x"]]]], + y = .data[[attributes(df)[["y"]]]] + ), + colour = "red", size = 1) + } + return(p) +} + +# This methods coors the backrounds based on the point density. This creates +# "landscape plots". +.add_point_density <- function(p, df, adjust = 1, ...){ + # Calculate the correct smoothing bandwidth. We could use the default + # values, but this custom bandwidth helps us to highlight the density + # better. + bandwidth <- adjust * c( + .get_bandwidth(df[[attributes(df)[["x"]]]]), + .get_bandwidth(df[[attributes(df)[["y"]]]])) + # Add the landscape + p <- p + + stat_density_2d( + aes(fill = after_stat(density)), + geom = "raster", + h = bandwidth, + contour = FALSE, + alpha = 0.5, + ) + + scale_fill_gradient(name = "Density", low = "white", high = "black") + return(p) +} + +# This function calculates the smoothing bandwidth. It highlights better point- +# rich areas than the default choice. +.get_bandwidth <- function(x){ + r <- quantile(x, c(0.25, 0.75)) + 4 * 1.06 * min(sd(x), (r[[2]] - r[[1]])/1.34) * length(x)^(-0.2) +} + +# This function adjusts the theme of the plot +.adjust_ordination_theme <- function( + p, df, + xlab = attributes(df)[["xlab"]], + ylab = attributes(df)[["ylab"]], + coord.equal = TRUE, + aspect.ratio = c("equal", "eigen", "free"), + ...){ + if( !.is_a_string(xlab) ){ + stop("'xlab' must be a single character value.", call. = FALSE) + } + if( !.is_a_string(ylab) ){ + stop("'ylab' must be a single character value.", call. = FALSE) + } + if( !.is_a_bool(coord.equal) ){ + stop("'coord.equal' must be TRUE or FALSE.", call. = FALSE) + } + aspect.ratio <- match.arg(aspect.ratio) + # + p <- p + theme_classic() + p <- p + labs(x = xlab, y = ylab) + if( !is.null(attributes(df)[["fill.by"]]) ){ + p <- p + labs(fill = attributes(df)[["fill.by"]]) + } + if( !is.null(attributes(df)[["colour.by"]]) ){ + p <- p + labs(colour = attributes(df)[["colour.by"]]) + } + if( !is.null(attributes(df)[["shape.by"]]) ){ + p <- p + labs(shape = attributes(df)[["shape.by"]]) + } + if( !is.null(attributes(df)[["size.by"]]) ){ + p <- p + labs(shape = attributes(df)[["size.by"]]) + } + if( !is.null(attributes(df)[["linetype.by"]]) ){ + p <- p + labs(linetype = attributes(df)[["linetype.by"]]) + } + + # Adjust axis length ratios + eigen_names <- c("eig", "percentVar") + if( aspect.ratio == "equal" ){ + # Use equal scaling on both axes so that one unit on the x-axis has the + # same physical length as one unit on the y-axis. This preserves + # distances and angles in the ordination and avoids visual distortion. + p <- p + coord_equal() + } else if ( aspect.ratio == "eigen" && + any(eigen_names %in% names(orig_attributes)) ){ + # Make the panel dimensions proportional to the eigenvalues so that the + # relative lengths of the axes reflect the variation explained by each + # ordination axis. + eigen_names <- eigen_names[ + eigen_names %in% names(orig_attributes)][[1L]] + eigen <- orig_attributes[[eigen_names]] + p <- p + theme(aspect.ratio = eigen[2] / eigen[1]) + } + + # Adjust colors + if( !is.null(attributes(df)[["colour.by"]]) ){ + name <- attributes(df)[["colour.by"]] + vals <- df[[name]] + p <- .resolve_plot_colours( + p, vals, name, + ) + } + if( !is.null(attributes(df)[["fill.by"]]) ){ + name <- attributes(df)[["fill.by"]] + vals <- df[[name]] + p <- .resolve_plot_colours( + p, vals, name, fill = TRUE + ) + } + + return(p) +} diff --git a/R/utils.R b/R/utils.R index 589e4429..429ebb9c 100644 --- a/R/utils.R +++ b/R/utils.R @@ -18,12 +18,15 @@ .is_function <- mia:::.is_function .get_name_in_parent <- mia:::.get_name_in_parent .is_an_integer <- mia:::.is_an_integer +.is_integer <- mia:::.is_integer TAXONOMY_RANKS <- mia:::TAXONOMY_RANKS .is_a_numeric <- mia:::.is_a_numeric .capitalize <- mia:::.capitalize .check_rowTree_present <- mia:::.check_rowTree_present .check_colTree_present <- mia:::.check_colTree_present .merge_features <- mia:::.merge_features +.check_dimred_present <- mia:::.check_dimred_present +.check_metadata_present <- mia:::.check_metadata_present .norm_label <- function(label, x){ if(!is.null(label)){ @@ -56,8 +59,9 @@ TAXONOMY_RANKS <- mia:::TAXONOMY_RANKS } # This function checks whether variable can be found from colData or rowData. +# Optionally, we can look also from rownames. .check_metadata_variable <- function( - tse, var, row = FALSE, col = FALSE, multiple = FALSE, + tse, var, row = FALSE, col = FALSE, multiple = FALSE, rownames = FALSE, var.name = .get_name_in_parent(var)){ if( !.is_a_bool(multiple) ){ stop("'multiple' must be TRUE or FALSE.", call. = FALSE) @@ -69,11 +73,14 @@ TAXONOMY_RANKS <- mia:::TAXONOMY_RANKS check_values <- c() check_values <- c(check_values, if(col) colnames(colData(tse))) check_values <- c(check_values, if(row) colnames(rowData(tse))) + msg_values <- check_values + # Optionally, we look also rownames + check_values <- c(check_values, if(rownames) rownames(tse)) var_found <- all( var %in% check_values ) if( !(is_string && var_found) ){ stop("'", var.name, "' must be ", ifelse(multiple, "", "a single "), "character value from the following options: '", - paste0(check_values, collapse = "', '"), "'", call. = FALSE) + paste0(msg_values, collapse = "', '"), "'", call. = FALSE) } } return(NULL) diff --git a/man/plotCCA.Rd b/man/plotCCA.Rd index 5a916a02..69572aa4 100644 --- a/man/plotCCA.Rd +++ b/man/plotCCA.Rd @@ -4,9 +4,7 @@ \alias{plotCCA} \alias{plotRDA} \alias{plotCCA,SingleCellExperiment-method} -\alias{plotCCA,matrix-method} \alias{plotRDA,SingleCellExperiment-method} -\alias{plotRDA,matrix-method} \title{Plot RDA or CCA object} \usage{ plotCCA(x, ...) @@ -15,17 +13,12 @@ plotRDA(x, ...) \S4method{plotCCA}{SingleCellExperiment}(x, dimred, ...) -\S4method{plotCCA}{matrix}(x, ...) - \S4method{plotRDA}{SingleCellExperiment}(x, dimred, ...) - -\S4method{plotRDA}{matrix}(x, ...) } \arguments{ \item{x}{a \code{\link[TreeSummarizedExperiment:TreeSummarizedExperiment-constructor]{TreeSummarizedExperiment}} -or a matrix of weights. The latter is returned as output from -\code{\link[mia:runCCA]{getRDA}}.} +object.} \item{...}{additional parameters for plotting, inherited from \code{\link[scater:plotReducedDim]{plotReducedDim}}, @@ -139,14 +132,17 @@ tse <- addRDA( formula = assay ~ ClinicalStatus + Gender + Age, distance = "bray", na.action = na.exclude - ) +) suppressWarnings({ # Create RDA plot coloured by variable plotRDA(tse, "RDA", colour.by = "ClinicalStatus") -# Create RDA plot with empty ellipses -plotRDA(tse, "RDA", colour.by = "ClinicalStatus", add.ellipse = "colour") +# Create RDA plot with ellipses +plotRDA( + tse, "RDA", colour.by = "ClinicalStatus", fill.by = "ClinicalStatus", + add.ellipse = TRUE +) # Create RDA plot with text encased in labels plotRDA(tse, "RDA", colour.by = "ClinicalStatus", vec.text = FALSE) @@ -156,18 +152,11 @@ plotRDA(tse, "RDA", colour.by = "ClinicalStatus", repel.labels = FALSE) # Create RDA plot without vectors plotRDA(tse, "RDA", colour.by = "ClinicalStatus", add.vectors = FALSE) - -# Calculate RDA as a separate object -rda_mat <- getRDA( - tse, - assay.type = "relabundance", - formula = assay ~ ClinicalStatus + Gender + Age, - distance = "bray", - na.action = na.exclude - ) - -# Create RDA plot from RDA matrix -plotRDA(rda_mat) }) } +\seealso{ +\itemize{ +\item \code{\link[=plotOrdination]{plotOrdination}} +} +} diff --git a/man/plotJointRPCA.Rd b/man/plotJointRPCA.Rd new file mode 100644 index 00000000..7a751124 --- /dev/null +++ b/man/plotJointRPCA.Rd @@ -0,0 +1,123 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/AllGenerics.R, R/plotJointRPCA.R +\name{plotJointRPCA} +\alias{plotJointRPCA} +\alias{plotJointRPCA,MultiAssayExperiment-method} +\alias{plotJointRPCA,SingleCellExperiment-method} +\title{Visualize Joint-RPCA results} +\usage{ +plotJointRPCA(x, ...) + +\S4method{plotJointRPCA}{MultiAssayExperiment}(x, dimred, ...) + +\S4method{plotJointRPCA}{SingleCellExperiment}(x, dimred, add.vectors = TRUE, ...) +} +\arguments{ +\item{x}{A +\code{\link[SummarizedExperiment:SummarizedExperiment-class]{SummarizedExperiment}} +or +\code{\link[MultiAssayExperiment:MultiAssayExperiment-class]{MultiAssayExperiment}} +object containing Joint-RPCA results.} + +\item{...}{Additional arguments passed to +\code{\link[=plotOrdination]{plotOrdination}} and to the feature-vector +plotting. +Commonly used Joint-RPCA-specific arguments include: +\itemize{ +\item \code{ntop}: \code{integer(1)} or \code{NULL}. Maximum number of +loading vectors shown per data layer. The vectors with the largest lengths +are retained. Set to \code{NULL} to display all vectors. (Default: \code{10}) +}} + +\item{dimred}{\code{character(1)}: Specifies the name of the Joint-RPCA +result to plot.} + +\item{add.vectors}{\code{logical(1)} or \code{character}. If +\code{TRUE}, feature loading vectors are added. Alternatively, a character +vector can be supplied to display only features whose names match the given +pattern(s). (Default: \code{TRUE})} +} +\value{ +A \code{ggplot2} object. +} +\description{ +Creates a two-dimensional ordination plot from Joint-RPCA results stored in a +\code{SingleCellExperiment} or \code{MultiAssayExperiment}. In addition to the +sample ordination, feature loadings can be visualized as vectors. +} +\details{ +This function is a wrapper around \code{\link[=plotOrdination]{plotOrdination}} +for Joint-RPCA results. Consequently, most graphical parameters are passed +directly to \code{plotOrdination()}, including options for colouring, +grouping, faceting and adding ellipses. See +\code{\link[=plotOrdination]{plotOrdination}} for a complete description of +these arguments. + +Feature loadings are plotted as vectors originating from the origin. By +default, only the longest loading vectors are shown for each data layer to +improve readability. +} +\examples{ +data("HintikkaXOData") +mae <- HintikkaXOData + +mae[[1]] <- transformAssay( + mae[[1]], + assay.type = "counts", + method = "rclr", + impute = FALSE +) +mae[[2]] <- transformAssay( + mae[[2]], + assay.type = "nmr", + method = "log10" +) + +mae <- addJointRPCA( + mae, + experiments = c(1, 2), + assay.types = c("rclr", "log10") +) + +# Basic Joint-RPCA plot +plotJointRPCA(mae, "JointRPCA") + +# Colour samples by metadata +plotJointRPCA(mae, "JointRPCA", colour.by = "Fat") + +# Add confidence ellipses +plotJointRPCA( + mae, + "JointRPCA", + colour.by = "Fat", + add.ellipse = TRUE +) + +# Show all loading vectors +plotJointRPCA( + mae, + "JointRPCA", + ntop = NULL +) + +# Display only selected feature vectors +plotJointRPCA( + mae, + "JointRPCA", + add.vectors = c("Bacteroides", "Roseburia") +) + +# Use boxed labels without repelling +plotJointRPCA( + mae, + "JointRPCA", + text.labels = FALSE, + repel.labels = FALSE +) + +} +\seealso{ +\itemize{ +\item \code{\link[=plotOrdination]{plotOrdination}} +} +} diff --git a/man/plotOrdination.Rd b/man/plotOrdination.Rd new file mode 100644 index 00000000..c005b8a9 --- /dev/null +++ b/man/plotOrdination.Rd @@ -0,0 +1,262 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/AllGenerics.R, R/plotOrdination.R +\name{plotOrdination} +\alias{plotOrdination} +\alias{plotOrdination,SingleCellExperiment-method} +\title{Visualize ordination results} +\usage{ +plotOrdination(x, ...) + +\S4method{plotOrdination}{SingleCellExperiment}(x, dimred, ...) +} +\arguments{ +\item{x}{a +\code{\link[SummarizedExperiment:SummarizedExperiment-class]{SummarizedExperiment}} +object.} + +\item{...}{Additional parameters for plotting. +\itemize{ +\item \code{ncomponents}: \code{integer vector} of length 2 or +\code{integer scalar}. Specifies which ordination components are plotted. +If a scalar is provided, the first two components are used. +(Default: \code{2L}) + +\item \code{colour.by}: \code{NULL} or \code{character scalar}. Specifies a +variable from \code{colData(x)} or a feature from \code{rownames(x)} used +to colour observations. Feature abundances are taken from +\code{assay.type}. (Default: \code{NULL}) + +\item \code{fill.by}: \code{NULL} or \code{character scalar}. Specifies a +variable from \code{colData(x)} or a feature from \code{rownames(x)} used +to fill observations or ellipses. Feature abundances are taken from +\code{assay.type}. Cannot be used together with +\code{add.density = TRUE}. (Default: \code{NULL}) + +\item \code{shape.by}: \code{NULL} or \code{character scalar}. Specifies a +categorical variable from \code{colData(x)} used for point shapes. +(Default: \code{NULL}) + +\item \code{size.by}: \code{NULL} or \code{character scalar}. Specifies a +variable from \code{colData(x)} used for point sizes. +(Default: \code{NULL}) + +\item \code{group.by}: \code{NULL} or \code{character scalar}. Specifies a +categorical variable from \code{colData(x)} used for grouping when drawing +ellipses, centroids and centroid vectors. (Default: \code{NULL}) + +\item \code{linetype.by}: \code{NULL} or \code{character scalar}. +Specifies a categorical variable from \code{colData(x)} used for ellipse +line types. (Default: \code{NULL}) + +\item \code{pair.by}: \code{NULL} or \code{character scalar}. Specifies a +variable from \code{colData(x)} identifying observations that should be +connected by lines. (Default: \code{NULL}) + +\item \code{sort.by}: \code{NULL} or \code{character scalar}. Specifies a +variable from \code{colData(x)} used to order observations before drawing +connecting lines. (Default: \code{NULL}) + +\item \code{facet.by}: \code{NULL} or \code{character scalar}. Specifies a +categorical variable from \code{colData(x)} used to split the plot into +facets. (Default: \code{NULL}) + +\item \code{assay.type}: \code{character scalar}. Name of the assay used +when \code{colour.by} or \code{fill.by} specifies a feature. +(Default: \code{"counts"}) + +\item \code{add.points}: \code{logical scalar}. Whether to draw sample +points. (Default: \code{TRUE}) + +\item \code{add.ellipse}: \code{logical scalar}. Whether to draw confidence +ellipses around groups. (Default: \code{FALSE}) + +\item \code{add.density}: \code{logical scalar}. Whether to draw a +two-dimensional density estimate in the background. +(Default: \code{FALSE}) + +\item \code{add.centroids}: \code{logical scalar}. Whether to draw group +centroids. (Default: \code{FALSE}) + +\item \code{add.centroids.lines}: \code{logical scalar}. Whether to connect +observations to their group centroids. (Default: \code{FALSE}) + +\item \code{add.vectors}: \code{logical scalar}. Whether to draw vectors +from the global centroid to group centroids. +(Default: \code{FALSE}) + +\item \code{add.rotation}: \code{logical scalar}. Whether to draw rotation +(species score) coordinates if available in the ordination result. +(Default: \code{add.species}) + +\item \code{add.species}: \code{logical scalar}. Alias for +\code{add.rotation}. (Default: \code{FALSE}) + +\item \code{add.expl.var}: \code{logical scalar}. Whether to append the +percentage of explained variance to the axis labels when available. +(Default: \code{FALSE}) + +\item \code{scales}: \code{character scalar}. Scaling used for faceted +plots. Passed to \code{ggplot2::facet_wrap()}. (Default: +\code{"fixed"}) + +\item \code{xlab}, \code{ylab}: \code{character scalar}. Axis labels. +Defaults to the ordination component names. + +\item \code{panel.by.eigen}: \code{logical scalar}. Whether to scale the +panel aspect ratio according to the eigenvalues of the ordination when +available. (Default: \code{TRUE}) + +\item \code{point.shape}: Shape used for points. +(Default: \code{19}) + +\item \code{point.alpha}: \code{numeric scalar}. Transparency of points. +Must be between 0 and 1. (Default: \code{0.4}) + +\item \code{ellipse.alpha}: \code{numeric scalar}. Transparency of ellipse +fills. Must be between 0 and 1. (Default: \code{0.2}) + +\item \code{ellipse.linewidth}: \code{numeric scalar}. Line width of +ellipse borders. (Default: \code{0.5} or \code{0} when +\code{fill.by} is specified.) + +\item \code{ellipse.linetype}: Integer specifying the ellipse line type. +(Default: \code{1}) + +\item \code{confidence.level}: \code{numeric scalar}. Confidence level used +for ellipse calculation. Must be between 0 and 1. +(Default: \code{0.95}) + +\item \code{adjust}: \code{numeric scalar}. Multiplicative adjustment for +the bandwidth used in the background density estimate. +(Default: \code{1}) +}} + +\item{dimred}{\code{character scalar}. Name of the reduced dimension result +stored in \code{reducedDim(x)} to visualize.} +} +\value{ +A \code{ggplot2} object. +} +\description{ +Creates a two-dimensional ordination plot from a reduced dimension result +stored in a +\code{\link[SingleCellExperiment:reducedDims]{SingleCellExperiment}} +object. Samples can be coloured, filled, shaped, sized, grouped or faceted +using sample metadata or feature abundances, and additional graphical +elements such as confidence ellipses, centroids, vectors and density +estimates can be added. +} +\details{ +This function provides a unified interface for visualizing ordination methods +such as PCA, PCoA, MDS, t-SNE, UMAP, RDA and CCA. The plotted coordinates are +retrieved from a reduced dimension result stored in +\code{reducedDim(x)}. +} +\examples{ +data("Tito2024QMP") +tse <- Tito2024QMP + +# Compute relative abundances and an MDS ordination +tse <- transformAssay(tse, method = "relabundance") +tse <- addMDS(tse, assay.type = "relabundance", method = "bray", ncomponents = 50) + +# Basic ordination plot +plotOrdination(tse, "MDS") + +# Colour samples by a sample-level variable +plotOrdination(tse, "MDS", colour.by = "diagnosis") + +# Colour samples by the abundance of a single feature +plotOrdination( + tse, "MDS", + colour.by = rownames(tse)[1], + assay.type = "relabundance") + +# Use multiple aesthetics simultaneously +plotOrdination( + tse, "MDS", + colour.by = "diagnosis", + shape.by = "colonoscopy" +) + +# Add confidence ellipses +plotOrdination( + tse, "MDS", + colour.by = "diagnosis", + group.by = "diagnosis", + add.ellipse = TRUE +) + +# Fill ellipses instead of colouring them +plotOrdination( + tse, "MDS", + fill.by = "diagnosis", + group.by = "diagnosis", + add.ellipse = TRUE +) + +# Show explained variance in the axis labels +plotOrdination( + tse, "MDS", + colour.by = "diagnosis", + add.expl.var = TRUE +) + +# Add group centroids +plotOrdination( + tse, "MDS", + colour.by = "diagnosis", + group.by = "diagnosis", + add.centroids = TRUE +) + +# Connect samples to their group centroids +plotOrdination( + tse, "MDS", + colour.by = "diagnosis", + group.by = "diagnosis", + add.centroids.lines = TRUE +) + +# Show vectors from the global centroid to each group centroid +plotOrdination( + tse, "MDS", + colour.by = "diagnosis", + group.by = "diagnosis", + add.vectors = TRUE +) + +# Draw a background density estimate +plotOrdination( + tse, "MDS", + add.density = TRUE +) + +# Split the plot into facets +plotOrdination( + tse, "MDS", + colour.by = "diagnosis", + facet.by = "colonoscopy" +) + +# Plot different ordination components +plotOrdination( + tse, "MDS", + ncomponents = c(2, 3) +) + +# Customize point appearance +plotOrdination( + tse, "MDS", + colour.by = "diagnosis", + point.shape = 17, + point.alpha = 0.8 +) + +} +\seealso{ +\itemize{ +\item \code{\link[scater:plotReducedDim]{scater::plotReducedDim}} +\item \code{\link[=plotCCA]{plotCCA}} +} +} diff --git a/tests/testthat/test-plotCCA.R b/tests/testthat/test-plotCCA.R index 798278ad..53393bdd 100644 --- a/tests/testthat/test-plotCCA.R +++ b/tests/testthat/test-plotCCA.R @@ -10,16 +10,15 @@ test_that("plot RDA/CCA", { # Run/calculate RDA tse <- addRDA(tse, assay.type = "counts", formula = assay ~ patient_status + cohort) - rda <- getRDA(tse, assay.type = "counts", formula = assay ~ patient_status + cohort) # Minimal functionality expect_no_error(plotRDA(tse, "RDA")) # Wrong-entry scenarios - expect_error(plotRDA(tse, "RDA", colour_by = "wrong colname")) - expect_error(plotRDA(tse, "RDA", colour_by = "cohort", shape_by = "wrong colname")) + expect_error(plotRDA(tse, "RDA", colour.by = "wrong colname")) + expect_error(plotRDA(tse, "RDA", colour.by = "cohort", shape.by = "wrong colname")) expect_error(plotRDA(tse, "RDA", add.ellipse = "invalid value"), - "'add.ellipse' must be one of c(TRUE, FALSE, 'fill', 'color').", + "'add.ellipse' must be TRUE or FALSE.", fixed = TRUE) expect_error(plotRDA(tse, "RDA", add.significance = "invalid value"), "'add.significance' must be TRUE or FALSE.") @@ -34,12 +33,12 @@ test_that("plot RDA/CCA", { ### 2). TEST plot layers ### - el_true <- plotRDA(tse, "RDA", colour_by = "patient_status") - el_false <- plotRDA(tse, "RDA", colour_by = "patient_status", add.ellipse = FALSE) - el_col <- plotRDA(tse, "RDA", colour_by = "patient_status", add.ellipse = "colour") - el_fill <- plotRDA(tse, "RDA", colour_by = "patient_status", add.ellipse = "fill") + el_true <- plotRDA(tse, "RDA", colour.by = "patient_status", add.ellipse = TRUE) + el_false <- plotRDA(tse, "RDA", colour.by = "patient_status", add.ellipse = FALSE) + el_col <- plotRDA(tse, "RDA", colour.by = "patient_status", add.ellipse = TRUE) + el_fill <- plotRDA(tse, "RDA", fill.by = "patient_status", add.ellipse = TRUE) expect_warning( - vec_false <- plotRDA(tse, "RDA", colour_by = "patient_status", add.vectors = FALSE) + vec_false <- plotRDA(tse, "RDA", colour.by = "patient_status", add.ellipse = TRUE, add.vectors = FALSE) ) # Filled ellipse has one more layer than no ellipse plot expect_equal(length(ggplot_build(el_true)[["data"]]), 4) @@ -51,7 +50,9 @@ test_that("plot RDA/CCA", { expect_false(all(ggplot_build(el_fill)[["data"]][[2]][["alpha"]] == 0)) # Check ggplot aesthetics - p_aes <- plotRDA(tse, "RDA", colour_by = "patient_status", ellipse.alpha = 0.5, + p_aes <- plotRDA(tse, "RDA", + colour.by = "patient_status", fill.by = "patient_status", + add.ellipse = TRUE, ellipse.alpha = 0.5, ellipse.linewidth = 0.2, ellipse.linetype = 3, vec.size = 0.6, vec.colour = "red", vec.linetype = 2, arrow.size = 0.15, label.colour = "blue", label.size = 5) @@ -72,8 +73,8 @@ test_that("plot RDA/CCA", { # expect_true(arrow_size == 0.15) # Vector or label text - p_vec <- plotRDA(tse, "RDA", colour_by = "patient_status", vec.text = TRUE) - p_lab <- plotRDA(tse, "RDA", colour_by = "patient_status", vec.text = FALSE) + p_vec <- plotRDA(tse, "RDA", colour.by = "patient_status", fill.by = "patient_status", add.ellipse = TRUE, vec.text = TRUE) + p_lab <- plotRDA(tse, "RDA", colour.by = "patient_status", fill.by = "patient_status", add.ellipse = TRUE, vec.text = FALSE) # There must be label column in both expect_true( "label" %in% (ggplot_build(p_vec)[["data"]][[4]] |> names()) ) expect_true( "label" %in% (ggplot_build(p_lab)[["data"]][[4]] |> names()) ) diff --git a/tests/testthat/test-plotJointRPCA.R b/tests/testthat/test-plotJointRPCA.R new file mode 100644 index 00000000..88100604 --- /dev/null +++ b/tests/testthat/test-plotJointRPCA.R @@ -0,0 +1,287 @@ + +.make_joint <- function() { + + tse <- makeTSE(nrow = 8, ncol = 4) + assayNames(tse) <- "counts" + + rd <- matrix(rnorm(8), ncol = 2) + class(rd) <- c("JointRPCA", class(rd)) + + rotation <- matrix(rnorm(16), ncol = 2) + rownames(rotation) <- paste0("feature", seq_len(nrow(rotation))) + + attr(rd, "rotation") <- rotation + attr(rd, "n_features") <- c(Layer1 = 4, Layer2 = 4) + + reducedDim(tse, "JointRPCA") <- rd + + return(tse) +} + +test_that("returns ggplot", { + + p <- plotJointRPCA(.make_joint(), "JointRPCA") + + expect_s3_class(p, "ggplot") + +}) + +test_that("FALSE disables vectors", { + + expect_no_error( + + plotJointRPCA( + .make_joint(), + "JointRPCA", + add.vectors = FALSE + ) + + ) + +}) + +test_that("character vectors accepted", { + + expect_no_error( + + plotJointRPCA( + .make_joint(), + "JointRPCA", + add.vectors = "feature" + ) + + ) + +}) + +test_that("invalid add.vectors errors", { + + expect_error( + + plotJointRPCA( + .make_joint(), + "JointRPCA", + add.vectors = 1 + ), + + "add.vectors" + + ) + +}) + +test_that("non JointRPCA reducedDim errors", { + + tse <- makeTSE(nrow = 5, ncol = 2) + + reducedDim(tse,"PCA") <- matrix(rnorm(4),2) + + expect_error( + + plotJointRPCA( + tse, + "PCA" + ), + + "JointRPCA" + + ) + +}) + +test_that("returns data.frame", { + + x <- .get_joint_rpca_vector_data( + .make_joint(), + "JointRPCA", + TRUE + ) + + expect_s3_class(x,"data.frame") + +}) + +test_that("FALSE returns NULL", { + + expect_null( + + .get_joint_rpca_vector_data( + .make_joint(), + "JointRPCA", + FALSE + ) + + ) + +}) + +test_that("feature filtering works", { + + x <- .get_joint_rpca_vector_data( + .make_joint(), + "JointRPCA", + "feature1" + ) + + expect_true( + + all(grepl( + "feature1", + x$vector_label + )) + + ) + +}) + +test_that("case insensitive filtering works", { + + x <- .get_joint_rpca_vector_data( + .make_joint(), + "JointRPCA", + "FEATURE", + ignore.case = TRUE + ) + + expect_gt(nrow(x),0) + +}) + +test_that("ntop limits each layer", { + + x <- .get_joint_rpca_vector_data( + .make_joint(), + "JointRPCA", + TRUE, + ntop = 2 + ) + + expect_true( + + all(table(x$Layer) <= 2) + + ) + +}) + +test_that("ntop NULL keeps all", { + + x <- .get_joint_rpca_vector_data( + .make_joint(), + "JointRPCA", + TRUE, + ntop = NULL + ) + + expect_equal( + + nrow(x), + + 8 + + ) + +}) + +test_that("ignore.case validated", { + + expect_error( + + .get_joint_rpca_vector_data( + .make_joint(), + "JointRPCA", + TRUE, + ignore.case = 1 + ), + + "ignore.case" + + ) + +}) + +test_that("ntop validated", { + + expect_error( + + .get_joint_rpca_vector_data( + .make_joint(), + "JointRPCA", + TRUE, + ntop = "a" + ), + + "ntop" + + ) + +}) + +.make_vectors <- function() { + + data.frame( + PC1 = c(.2,.5), + PC2 = c(.4,.1), + Layer = c("A","B"), + vector_label = c("x","y") + ) + +} + +test_that("returns ggplot", { + + p <- ggplot() + + p2 <- .add_ordination_vectors( + p, + .make_vectors() + ) + + expect_s3_class(p2,"ggplot") + + expect_length( + + p2$layers, + + 2 + + ) + +}) + +test_that("supports geom_label", { + + p <- .add_ordination_vectors( + ggplot(), + .make_vectors(), + text.labels = FALSE + ) + + expect_s3_class(p,"ggplot") + +}) + +test_that("supports non-repel labels", { + + p <- .add_ordination_vectors( + ggplot(), + .make_vectors(), + repel.labels = FALSE + ) + + expect_s3_class(p,"ggplot") + +}) + +test_that("supports label without repel", { + + p <- .add_ordination_vectors( + ggplot(), + .make_vectors(), + repel.labels = FALSE, + text.labels = FALSE + ) + + expect_s3_class(p,"ggplot") + +}) diff --git a/tests/testthat/test-plotOrdination.R b/tests/testthat/test-plotOrdination.R new file mode 100644 index 00000000..4b85796a --- /dev/null +++ b/tests/testthat/test-plotOrdination.R @@ -0,0 +1,228 @@ + +.make_tse <- function() { + tse <- makeTSE(ncol = 2, nrow = 10) + assayNames(tse) <- "counts" + + reducedDim(tse, "PCA") <- matrix( + seq_len(6), + ncol = 3, + dimnames = list(colnames(tse), c("PC1", "PC2", "PC3")) + ) + + attr(reducedDim(tse, "PCA"), "eig") <- c(60, 30, 10) + + colData(tse)$numeric <- c(1, 2) + colData(tse)$time <- c(1, 2) + + return(tse) +} + +test_that("defaults are returned", { + + tse <- .make_tse() + + args <- .check_ordination_input(tse, "PCA") + + expect_identical(args$x, tse) + expect_identical(args$dimred, "PCA") + expect_identical(args$ncomponents, 1:2) + expect_true(args$add.expl.var == FALSE) + +}) + +test_that("scalar ncomponents expands", { + + args <- .check_ordination_input( + .make_tse(), + "PCA", + ncomponents = 2L + ) + + expect_identical(args$ncomponents, c(1L, 2L)) + +}) + +test_that("vector ncomponents is preserved", { + + args <- .check_ordination_input( + .make_tse(), + "PCA", + ncomponents = c(2L, 3L) + ) + + expect_identical(args$ncomponents, c(2L, 3L)) + +}) + +test_that("three components are rejected", { + + expect_error( + .check_ordination_input( + .make_tse(), + "PCA", + ncomponents = 1:3 + ), + "ncomponents" + ) + +}) + +test_that("component outside range is rejected", { + + expect_error( + .check_ordination_input( + .make_tse(), + "PCA", + ncomponents = c(1L, 4L) + ), + "ncomponents" + ) + +}) + +test_that("non-integer components are rejected", { + + expect_error( + .check_ordination_input( + .make_tse(), + "PCA", + ncomponents = c(1, 2.5) + ), + "ncomponents" + ) + +}) + +test_that("density and fill are mutually exclusive", { + + expect_error( + .check_ordination_input( + .make_tse(), + "PCA", + fill.by = "group", + add.density = TRUE + ), + "Both 'add.density' and 'fill.by' cannot be specified simultaneously." + ) + +}) + +for (arg in c( + "add.points", + "add.ellipse", + "add.density", + "add.centroids", + "add.centroids.lines", + "add.vectors", + "add.rotation", + "add.expl.var" +)) { + + test_that(paste(arg, "must be logical"), { + + x <- list( + x = .make_tse(), + dimred = "PCA" + ) + + x[[arg]] <- 1 + + expect_error( + do.call(.check_ordination_input, x), + arg + ) + + }) + +} + +test_that("extra arguments are propagated", { + + args <- .check_ordination_input( + .make_tse(), + "PCA", + point.alpha = .5, + point.shape = 17 + ) + + expect_identical(args$point.alpha, .5) + expect_identical(args$point.shape, 17) + +}) + +test_that("metadata arguments are propagated", { + + args <- .check_ordination_input( + .make_tse(), + "PCA", + colour.by = "group", + fill.by = "group", + shape.by = "group", + size.by = "numeric", + group.by = "group", + linetype.by = "group", + pair.by = "ID", + sort.by = "time", + facet.by = "group" + ) + + expect_identical(args$colour.by, "group") + expect_identical(args$fill.by, "group") + expect_identical(args$shape.by, "group") + expect_identical(args$size.by, "numeric") + expect_identical(args$group.by, "group") + expect_identical(args$linetype.by, "group") + expect_identical(args$pair.by, "ID") + expect_identical(args$sort.by, "time") + expect_identical(args$facet.by, "group") + +}) + +test_that("plotOrdination returns a ggplot", { + p <- plotOrdination(.make_tse(), "PCA") + expect_s3_class(p, "ggplot") +}) + +test_that("integer ncomponents of length one expands to first two components", { + + args <- .check_ordination_input( + .make_tse(), + "PCA", + ncomponents = 2L + ) + + expect_identical(args$ncomponents, c(1L, 2L)) + +}) + +test_that("component order is preserved", { + + args <- .check_ordination_input( + .make_tse(), + "PCA", + ncomponents = c(3L, 1L) + ) + + expect_identical(args$ncomponents, c(3L, 1L)) + +}) + +test_that("all plotting flags can be enabled", { + + expect_no_error( + + .check_ordination_input( + .make_tse(), + "PCA", + add.points = TRUE, + add.ellipse = TRUE, + add.centroids = TRUE, + add.centroids.lines = TRUE, + add.vectors = TRUE, + add.rotation = TRUE, + add.expl.var = TRUE + ) + + ) + +})