Add a function for using PCA to summarize expression across several g…

…enes before plotting spatially

NAMESPACE

@@ -1,6 +1,7 @@
 # Generated by roxygen2: do not edit by hand
 
 export(spot_plot)
+export(spot_plot_pca)
 export(spot_plot_sparsity)
 export(spot_plot_z_score)
 import(Matrix)

R/multi_gene.R

@@ -124,6 +124,78 @@ spot_plot_sparsity <- function(spe, genes, sample_id, assayname = "counts", minC
     return(p)
 }
 
+#' Spatial plot of multiple genes combined by averaging Z scores
+#'
+#' This function wraps around \code{spot_plot}, plotting a summary of expression
+#' across multiple genes in a single spatial plot. The specified expression
+#' assay is subsetted to the select genes and transposed before performing PCA.
+#' Each PC thus represents a "spot profile", where larger elements represent
+#' spots of greater expression variation within the gene set. The first such
+#' PC is plotted spatially, negated if necessary so that the majority of
+#' coefficients across genes are positive (this encourages higher values to
+#' represent areas of higher expression of the genes).
+#'
+#' @inheritParams spot_plot_z_score
+#'
+#' @return A \code{ggplot} object containing a "spot plot" of the specified
+#' sample and genes
+#'
+#' @export
+#' @author Nicholas J. Eagles
+#' @import SpatialExperiment SummarizedExperiment Matrix
+#' @family Spot plots summarizing expression of multiple genes simultaneously
+#'
+#' @examples
+#'
+#' #   Grab an example SpatialExperiment and suppose all of its spots should be
+#' #   plotted (for spatialNAc, 'exclude_overlapping' will only have genuinely
+#' #   overlapping spots be TRUE)
+#' spe <- spatialLIBD::fetch_data(type = "spatialDLPFC_Visium_example_subset")
+#' spe$exclude_overlapping <- FALSE
+#'
+#' white_matter_genes <- c(
+#'     "ENSG00000197971", "ENSG00000131095", "ENSG00000123560",
+#'     "ENSG00000171885"
+#' )
+#' spot_plot_pca(
+#'     spe = spe,
+#'     genes = white_matter_genes,
+#'     sample_id = unique(spe)$sample_id[1]
+#' )
+spot_plot_pca <- function(spe, genes, sample_id, assayname = "logcounts", minCount = 0, ...) {
+    #   Check validity of arguments
+    .multi_gene_validity_check(
+        spe, genes, sample_id, assayname, minCount, ...
+    )
+
+    spe <- spe[genes, spe$sample_id == sample_id]
+
+    pc_exp = prcomp(t(assays(spe)[[assayname]]), center = TRUE, scale = TRUE)
+    spe$pc_select_genes <- pc_exp$x[,'PC1']
+
+    #   Given that:
+    #       - 'genes' is assumed to represent markers of the subregion (and
+    #         thus their expression follows a similar pattern spatially)
+    #       - the first PC captures the most variation
+    #   Then each gene's coefficients to the first PC should tend to have
+    #   the same sign. Next, the sign of each PC is arbitary, and we'd like
+    #   plots to have positive values where expression is greater. If most
+    #   genes have negative coefficients to the first PC, we reverse the
+    #   sign of the coefficients to make visual intrepretation consistent
+    if (mean(pc_exp$rotation[,1] > 0) < 0.5) {
+        spe$pc_select_genes = -1 * spe$pc_select_genes
+    }
+
+    #   Plot spatial distribution of this proportion
+    p <- spot_plot(
+        spe, sample_id,
+        var_name = "pc_select_genes",
+        is_discrete = FALSE, minCount = minCount, assayname = assayname, ...
+    )
+
+    return(p)
+}
+
 #'   Check the validity of arguments passed to \code{multi_gene.R} plotting functions
 #'
 #' @author Nicholas J. Eagles

R/spot_plot.R

@@ -38,30 +38,29 @@
 #' @export
 #' @author Nicholas J. Eagles
 #' @import viridisLite spatialLIBD ggplot2 SpatialExperiment SummarizedExperiment
-#' 
-#' @examples 
-#' 
+#'
+#' @examples
+#'
 #' #   Grab an example SpatialExperiment and suppose all of its spots should be
 #' #   plotted (for spatialNAc, 'exclude_overlapping' will only have genuinely
 #' #   overlapping spots be TRUE)
 #' spe <- spatialLIBD::fetch_data(type = "spatialDLPFC_Visium_example_subset")
 #' spe$exclude_overlapping <- FALSE
 #'
 #' #   Plot age spatially for the first sample
-#' sample_id = unique(spe$sample_id)[1]
+#' sample_id <- unique(spe$sample_id)[1]
 #' p <- spot_plot(
-#'    spe, sample_id = sample_id,
-#'    title = sample_id, var_name = "age",
-#'    include_legend = TRUE, is_discrete = FALSE, minCount = 0,
-#'    assayname = "logcounts"
+#'     spe,
+#'     sample_id = sample_id,
+#'     title = sample_id, var_name = "age",
+#'     include_legend = TRUE, is_discrete = FALSE, minCount = 0,
+#'     assayname = "logcounts"
 #' )
 #' print(p)
-spot_plot <- function(
-        spe, sample_id, image_id = "lowres",
-        title = sprintf('%s_%s', sample_id, var_name), var_name,
-        include_legend = TRUE, is_discrete, colors = NULL,
-        assayname = "logcounts", minCount = 0.5
-    ) {
+spot_plot <- function(spe, sample_id, image_id = "lowres",
+    title = sprintf("%s_%s", sample_id, var_name), var_name,
+    include_legend = TRUE, is_discrete, colors = NULL,
+    assayname = "logcounts", minCount = 0.5) {
     #   This value was determined empirically, and results in good spot sizes.
     #   Note that it's sample-independent, and the final spot size to pass to
     #   'vis_gene' or 'vis_clus' uses this value along with the image

tests/testthat/test-multi_gene.R

@@ -54,3 +54,6 @@ test_that("spot_plot_z_score", {
 test_that("spot_plot_sparsity", {
     test_plot_fun(spe, spot_plot_sparsity)
 })
+test_that("spot_plot_pca", {
+    test_plot_fun(spe, spot_plot_pca)
+})