Fix ANOVA and pairwise gene set enrichment and spatial registration plot visualization

DESCRIPTION

@@ -80,7 +80,8 @@ Imports:
     dplyr,
     ComplexHeatmap,
     circlize,
-    scrapper
+    scrapper,
+    scales
 Roxygen: list(markdown = TRUE)
 URL: https://github.com/LieberInstitute/spatialLIBD
 BugReports: https://support.bioconductor.org/tag/spatialLIBD

R/app_server.R

@@ -1540,7 +1540,16 @@ app_server <- function(input, output, session) {
             assayname = input$layer_model_assayname
         )
     })
-
+
+    # handle anova single colum issue
+    modeling_results_app <- modeling_results
+
+    # check anova col names, those ending with '_all" unchanged
+    if ("anova" %in% names(modeling_results_app)) {
+        colnames(modeling_results_app$anova) <-
+            sub("_$", "_all", colnames(modeling_results_app$anova))
+    }
+
     static_layer_gene_set_enrichment <- reactive({
         if (!is.null(input$geneSet)) {
             gene_list <-
@@ -1572,12 +1581,13 @@ app_server <- function(input, output, session) {
             )
             gene_list <- asd_sfari_geneList
         }
-
+
+        ## single-column ANOVA results being dropped to vectors
         ## Run the enrichment
         gene_set_enrichment(
             gene_list,
             input$layer_gene_fdrcut,
-            modeling_results,
+            modeling_results_app,
             input$layer_model,
             reverse = ifelse(
                 input$layer_model != "anova",
@@ -1589,18 +1599,54 @@ app_server <- function(input, output, session) {
 
     static_layer_gene_set_enrichment_plot <- reactive({
         enriched <- static_layer_gene_set_enrichment()
+
+        # Convert SetSize from a list column to numeric
+        enriched$SetSize <- as.numeric(unlist(
+            enriched$SetSize, use.names = FALSE)
+            )
+
+        # Original domain colors
+        model_colors <- sce_layer[[paste0(default_cluster, "_colors")]]
+
+        # if no default_cluster colors- use colors from spatialLIBD::get_colors()
+        if (is.null(model_colors)) {
+            model_colors <- spatialLIBD::get_colors(
+                clusters = unique(enriched$test),
+                use_libd_layer_colors = TRUE
+                )
+        }
+
+        # Add colors for pairwise tests
+        pairwise_tests <- unique(
+            enriched$test[enriched$model_type == "pairwise"]
+            )
+
+        for (test in pairwise_tests) {
+            if (!test %in% names(model_colors)) {
+                groups <- strsplit(test, "-", fixed = TRUE)[[1]]
+
+                if (all(groups %in% names(model_colors))) {
+                    model_colors[test] <- scales::colour_ramp(
+                        model_colors[groups]
+                        )(0.5) #midpoint between 2 input colors
+                } else {
+                    model_colors[test] <- "grey60"
+                }
+            }
+        }
+
+        # ANOVA has one overall test, not a spatial-domain-specific test
+        if (input$layer_model == "anova") {
+            model_colors <- "grey60" 
+            names(model_colors) <- unique(as.character(enriched$test))
+
+        }
         gene_set_enrichment_plot(
             enriched,
             plot_SetSize_bar = TRUE,
-            model_colors = get_colors(
-                if (!is.null(sce_layer[[paste0(default_cluster, "_colors")]])) {
-                    sce_layer[[paste0(default_cluster, "_colors")]]
-                } else {
-                    spatialLIBD::libd_layer_colors
-                },
-                clusters = unique(enriched$test)
-            )
+            model_colors = model_colors[unique(enriched$test)]
         )
+
     })
 
     static_layer_external_tstat <- reactive({
@@ -1620,8 +1666,24 @@ app_server <- function(input, output, session) {
                 spatialLIBD::tstats_Human_DLPFC_snRNAseq_Nguyen_topLayer
         }
 
+        cor_results <- modeling_results_app
+
+        ## single-column ANOVA results being dropped to vectors
+        if (input$layer_model == "anova") {
+            cor_results$anova$f_stat_tmp <-
+                cor_results$anova$f_stat_all
+        }
+
         ## Compute the correlations
-        layer_stat_cor(input_stat, modeling_results, input$layer_model)
+        cor_mat <- layer_stat_cor(input_stat, cor_results, input$layer_model)
+
+        ## Remove the temporary duplicate column from the output
+        if (input$layer_model == "anova") {
+            cor_mat <- cor_mat[, colnames(cor_mat) != "tmp", drop = FALSE]
+        }
+
+        cor_mat
+
     })
 
     static_layer_external_tstat_annotated_clusters <- reactive({
@@ -1633,26 +1695,53 @@ app_server <- function(input, output, session) {
     })
 
     static_layer_external_tstat_plot <- reactive({
+        cor_mat <- static_layer_external_tstat()
+
+        # Get colors for individual spatial domains
+        model_colors <- sce_layer[[paste0(default_cluster, "_colors")]]
+
+        if (is.null(model_colors)) {
+            model_colors <- spatialLIBD::get_colors(
+                clusters = colnames(cor_mat),
+                use_libd_layer_colors = TRUE
+                )
+        }
+
+        # Add colors for pairwise comparisons, like L1-L2
+        pairwise_tests <- colnames(cor_mat)[grepl("-", colnames(cor_mat), fixed = TRUE)]
+
+        for (test in pairwise_tests) {
+            if (!test %in% names(model_colors)) {
+                groups <- strsplit(test, "-", fixed = TRUE)[[1]]
+
+                if (all(groups %in% names(model_colors))) {
+                    model_colors[test] <- scales::colour_ramp(
+                        model_colors[groups]
+                        )(0.5) ## midpoint between the two input colors
+                } else {
+                    model_colors[test] <- "grey60"
+                }
+            }
+        }
+
+        # ANOVA has one overall test, not a spatial-domain-specific test
+        if (input$layer_model == "anova") {
+            model_colors <- "grey60"
+            names(model_colors) <- colnames(cor_mat)
+        }
+
         layer_stat_cor_plot(
-            static_layer_external_tstat(),
+            cor_mat,
             max(c(0.1, input$layer_tstat_max)),
             -max(c(0.1, input$layer_tstat_max)),
             cluster_rows = FALSE,
             cluster_columns = FALSE,
             annotation = static_layer_external_tstat_annotated_clusters(),
-            query_colors = get_colors(
-                clusters = rownames(static_layer_external_tstat())
-            ),
-            reference_colors = get_colors(
-                if (!is.null(sce_layer[[paste0(default_cluster, "_colors")]])) {
-                    sce_layer[[paste0(default_cluster, "_colors")]]
-                } else {
-                    spatialLIBD::libd_layer_colors
-                },
-                clusters = colnames(static_layer_external_tstat())
-            )
+            query_colors = get_colors(clusters = rownames(cor_mat)),
+            reference_colors = model_colors
         )
     })
+
 
     ## layer download PDF buttons
     output$layer_downloadReducedDim <- downloadHandler(

R/gene_set_enrichment.R

@@ -101,7 +101,7 @@ gene_set_enrichment <-
         }
 
         tstats <-
-            model_results[, grep("[f|t]_stat_", colnames(model_results))]
+            model_results[, grep("[f|t]_stat_", colnames(model_results)), drop=FALSE]
         colnames(tstats) <-
             gsub("[f|t]_stat_", "", colnames(tstats))
 
@@ -125,7 +125,7 @@ gene_set_enrichment <-
         }
 
         fdrs <-
-            model_results[, grep("fdr_", colnames(model_results))]
+            model_results[, grep("fdr_", colnames(model_results)), drop=FALSE]
 
         enrichTab <-
             do.call(

R/gene_set_enrichment_plot.R

@@ -212,7 +212,7 @@ gene_set_enrichment_plot <-
             colnames(res) <-
                 gsub(".*_mypattern_", "", colnames(res))
             rownames(res) <- unique(enrichment$ID)
-            res <- res[, levels(as.factor(enrichment$test))]
+            res <- res[, levels(as.factor(enrichment$test)), drop=FALSE]
             t(res)
         }
         wide_or <- make_wide("OR_char")

R/get_colors.R

@@ -9,6 +9,8 @@
 #' generate up to 36 unique colors. If the number of unique clusters is beyond
 #' 36 then this function will fail.
 #' @param clusters A vector of cluster names.
+#' @param use_libd_layer_colors Logical. If `TRUE` and `colors` is `NULL`,
+#' use `libd_layer_colors` when all cluster names match its names.
 #'
 #' @return A named vector where the values are the colors to use for
 #' displaying them different clusters. For some use cases, you might have to
@@ -41,9 +43,19 @@
 #' )
 #' log_var_sorted <- sort_clusters(log_var)
 #' get_colors(colors = NULL, clusters = log_var_sorted)
-get_colors <- function(colors = NULL, clusters) {
+get_colors <- function(colors = NULL, clusters, use_libd_layer_colors = FALSE) {
     n_clus <- length(unique(clusters[!is.na(clusters)]))
 
+    # Use predefined spatialLIBD layer colors when requested
+    if (
+        use_libd_layer_colors &&
+        is.null(colors) &&
+        all(unique(clusters) %in% names(libd_layer_colors))
+    ) {
+        colors <- libd_layer_colors
+    }
+
+    ## if colors is NULL or user supplies insufficent colors
     if (is.null(colors) | n_clus > length(colors)) {
         ## Original ones
         # colors <- c("#b2df8a","#e41a1c","#377eb8","#4daf4a","#ff7f00","gold",

R/run_app.R

@@ -10,6 +10,23 @@
 #' will disable the pseudo-bulked portion of the web application. See the
 #' examples for one such case as well as the vignette that describes how
 #' you can use `spatialLIBD` with public data sets provided by 10x Genomics.
+#' 
+#'  ---- Optional: Add custom cluster colors ----
+#'
+#' Custom colors can be specified by adding a color column to the
+#' `colData` of the input objects. The column name must be the cluster
+#' variable name followed by `_colors` (for example, `BayesSpace_colors`
+#' when `default_cluster = "BayesSpace"`).
+#'
+#' Define a named vector mapping cluster labels to colors, then assign
+#' the corresponding colors to each observation. Retain the cluster
+#' labels as names in the color vector.
+#'
+#' When using layer-level data (`sce_layer`), the same color column
+#' can be added to that object. Colors should correspond to the
+#' cluster labels in the selected `default_cluster` variable.
+#' If no custom colors are provided, the app uses `spatialLIBD::get_colors()`.
+#' 
 #' That vignette is available at
 #' <http://research.libd.org/spatialLIBD/articles/TenX_data_download.html>.
 #'
@@ -323,7 +340,7 @@ run_app <- function(
         )))
     }
 
-    with_golem_options(
+    golem::with_golem_options(
         app = shinyApp(ui = app_ui, server = app_server),
         golem_opts = list(
             spe = spe,

vignettes/spatialLIBD.Rmd

@@ -263,6 +263,54 @@ dim(sig_genes)
 
 Now that you have downloaded the data, you can interactively explore the data using a `r CRANpkg('shiny')` `r Citep(bib[['shiny']])` web application contained within `r Biocpkg('spatialLIBD')` `r Citep(bib[['spatialLIBD']])`. To do so, use the `run_app()` function as shown below:
 
+### Customizing cluster colors
+
+The `spatialLIBD` application supports custom colors for spatial domains or clusters in both the spot-level (`spe`) and layer-level (`sce_layer`) data.
+
+To specify custom colors, add a color column to the `colData()` of the corresponding input object. The column name must be the name of the cluster variable followed by `_colors`. For example, if the cluster variable is `layer_guess`, the corresponding color column must be named `layer_guess_colors`.
+
+Each color column contains one color value per observation in the corresponding object:
+- In `spe`, the colors are assigned to individual spots according to their cluster labels.
+- In `sce_layer`, the colors are assigned to the layer-level groups according to their cluster labels.
+
+The colors can be specified using a named vector, where the names correspond to the cluster labels and the values are valid R colors.
+
+For example, the following code assigns custom colors to the clusters in both objects:
+
+```{r custom_cluster_colors, eval = FALSE}
+## Define a named vector of colors for each cluster
+my_cluster_colors <- c(
+    Cluster1 = "#4E79A7",
+    Cluster2 = "#F28E2B",
+    Cluster3 = "#E15759",
+    Cluster4 = "#76B7B2",
+    Cluster5 = "#59A14F",
+    Cluster6 = "#EDC948",
+    Cluster7 = "#B07AA1",
+    Cluster8 = "#FF9DA7",
+    Cluster9 = "#9C755F"
+)
+
+## Specify the cluster variable used by the application
+cluster_var <- "my_cluster"
+
+## Add the corresponding color column to the spot-level object
+colData(spe)[[paste0(cluster_var, "_colors")]] <-
+    my_cluster_colors[
+        as.character(colData(spe)[[cluster_var]])
+    ]
+
+## Add the corresponding color column to the layer-level object
+colData(sce_layer)[[paste0(cluster_var, "_colors")]] <-
+    my_cluster_colors[
+        as.character(colData(sce_layer)[[cluster_var]])
+    ]
+```
+
+Custom color columns are optional. If they are not provided, the application uses the predefined `spatialLIBD::libd_layer_colors` palette when the cluster labels match its names. This palette contains colors for `Layer1` through `Layer6`, `WM`, `NA`, and `WM2`. For cluster labels that do not match these predefined names, the application generates a palette using `spatialLIBD::get_colors()`, with one color assigned to each cluster.
+
+The application uses the corresponding colors in visualizations of the spot-level and layer-level data, including spatial plots, reduced-dimension plots, gene set enrichment results and spatial registration plots for ANOVA, enrichment and pairwise model results.
+
 ```{r}
 if (interactive()) {
     run_app(