Add layer_stat_cor_plot_complex()

DESCRIPTION

@@ -79,7 +79,8 @@ Imports:
     statmod,
     MatrixGenerics,
     rlang,
-    dplyr
+    dplyr,
+    ComplexHeatmap
 RoxygenNote: 7.3.2
 Roxygen: list(markdown = TRUE)
 URL: https://github.com/LieberInstitute/spatialLIBD

NAMESPACE

@@ -56,6 +56,9 @@ import(plotly, except = last_plot)
 import(shiny)
 importFrom(AnnotationHub,query)
 importFrom(BiocGenerics,which)
+importFrom(ComplexHeatmap,Heatmap)
+importFrom(ComplexHeatmap,columnAnnotation)
+importFrom(ComplexHeatmap,rowAnnotation)
 importFrom(DT,DTOutput)
 importFrom(DT,renderDT)
 importFrom(GenomicRanges,seqnames)

R/annotate_registered_clusters.R

@@ -59,13 +59,6 @@ annotate_registered_clusters <-
                 cutoff_merge_ratio = cutoff_merge_ratio
             )
 
-        if (all(colnames(cor_stats_layer) %in% c("WM", paste0("Layer", seq_len(6))))) {
-            ## Simplify names when working with the default data
-            annotated <- gsub("ayer", "", annotated)
-            annotated <- gsub("\\/L", "\\/", annotated)
-            annotated <- gsub("^WM\\/", "WM\\/L", annotated)
-        }
-
         confidence <- apply(cor_stats_layer, 1, max) > confidence_threshold
 
         result <- data.frame(
@@ -83,6 +76,16 @@ annotate_registered_clusters <-
                 result$layer_label,
                 ifelse(result$layer_confidence == "good", "", "*")
             )
+
+        ## Add simplified label for WM/Layer annotations
+        if (all(colnames(cor_stats_layer) %in% c("WM", paste0("Layer", seq_len(6))))) {
+            result$layer_label_simple <- result$layer_label
+            ## Simplify names when working with the default data
+            result$layer_label_simple  <- gsub("ayer", "", result$layer_label_simple )
+            result$layer_label_simple  <- gsub("\\/L", "\\/", result$layer_label_simple )
+            result$layer_label_simple  <- gsub("^WM\\/", "WM\\/L", result$layer_label_simple )
+        }
+
         return(result)
     }
 

R/layer_stat_cor_plot.R

@@ -1,126 +1,182 @@
-#' Visualize the layer modeling correlation of statistics
-#'
-#' This function makes a heatmap from the [layer_stat_cor()] correlation matrix
-#' between a given set of cell cluster/type statistics derived from scRNA-seq
-#' or snRNA-seq data (among other types) and the layer statistics from the
-#' Human DLPFC Visium data (when using the default arguments).
-#'
+#' Visualize the correlation of layer modeling t-statistics with ComplexHeatmap
+#' 
+#' Use ComplexHeatmap to plot the correlation matrix between a reference and 
+#' query modeling statistics from [layer_stat_cor()]. 
+#' 
+#' Includes functionality to add color annotations, 
+#' (helpful to match to colors in Visium spot plots), and annotations from 
+#' [annotate_registered_clusters()].
+#' 
 #' @param cor_stats_layer The output of [layer_stat_cor()].
-#' @param max A `numeric(1)` specifying the highest correlation value for the
-#' color scale (should be between 0 and 1).
-#' @param min A `numeric(1)` specifying the lowest correlation value for the
-#' color scale (should be between 0 and -1).
-#' @param layerHeights A `numeric()` vector of length equal to
-#' `ncol(cor_stats_layer) + 1` that starts at 0 specifying where
-#' to plot the y-axis breaks which can be used for re-creating the length of
-#' each brain layer. Gets passed to [layer_matrix_plot()].
-#' @param cex Passed to [layer_matrix_plot()].
+#' @param color_max A `numeric(1)` specifying the highest correlation value for 
+#' the color scale (should be between 0 and 1).
+#' @param color_min A `numeric(1)` specifying the lowest correlation value for 
+#' the color scale (should be between 0 and -1).
+#' @param color_scale A `character` vector specifying the color scale for the 
+#' fill of the heatmap, defaults to classic purple -> green.
+#' @param query_colors named `character` vector of colors, Adds colors to query 
+#' row annotations.
+#' @param reference_colors named `character` vector of colors, Adds colors to 
+#' reference column annotations.
+#' @param annotation annotation data.frame output of [annotate_registered_clusters()],
+#' adds 'X' for good confidence annotations, '*' for poor confidence.
+#' @param ... Additional parameters passed to [ComplexHeatmap::Heatmap()][ComplexHeatmap::Heatmap()]
+#' ex. `cluster_rows` and `cluster_columns`.
 #'
-#' @return A heatmap for the correlation matrix between statistics.
+#'
+#' @return ([Heatmap-class][ComplexHeatmap::Heatmap-class]) plot of t-stat correlations
 #' @export
-#' @author Andrew E Jaffe, Leonardo Collado-Torres
+#' @author Louise Huuki-Myers
 #' @family Layer correlation functions
-#' @seealso layer_matrix_plot annotate_registered_clusters
 #'
 #' @importFrom RColorBrewer brewer.pal
 #' @importFrom grDevices colorRampPalette
-#' @details Check
-#' https://github.com/LieberInstitute/HumanPilot/blob/master/Analysis/Layer_Guesses/dlpfc_snRNAseq_annotation.R
-#' for a full analysis from which this family of functions is derived from.
+#' @importFrom ComplexHeatmap columnAnnotation rowAnnotation Heatmap
 #'
 #' @examples
-#'
 #' ## Obtain the necessary data
+#' ## reference human pilot modeling results
 #' if (!exists("modeling_results")) {
 #'     modeling_results <- fetch_data(type = "modeling_results")
 #' }
+#' 
+#' ## querey spatailDLPFC modeling 
+#' query_modeling_results <- fetch_data(type = "spatialDLPFC_Visium_modeling_results")
+#' 
+#' ## extract t-statics and rename
+#' registration_t_stats <- query_modeling_results$enrichment[, grep("^t_stat", colnames(query_modeling_results$enrichment))]
+#' colnames(registration_t_stats) <- gsub("^t_stat_", "", colnames(registration_t_stats))
 #'
 #' ## Compute the correlations
 #' cor_stats_layer <- layer_stat_cor(
-#'     tstats_Human_DLPFC_snRNAseq_Nguyen_topLayer,
+#'     stats = registration_t_stats,
 #'     modeling_results,
 #'     model_type = "enrichment"
 #' )
 #'
 #' ## Visualize the correlation matrix
-#' layer_stat_cor_plot(cor_stats_layer, max = max(cor_stats_layer))
-#'
-#' ## Annotate then re-plot
-#' rownames(cor_stats_layer) <- paste0(
-#'     rownames(cor_stats_layer),
-#'     " - ",
-#'     annotate_registered_clusters(cor_stats_layer)$layer_label
-#' )
-#' layer_stat_cor_plot(cor_stats_layer, max = max(cor_stats_layer))
-#'
-#' ## Restrict the range of colors further
-#' layer_stat_cor_plot(cor_stats_layer, max = 0.25)
+#' 
+#' ## most basic
+#' layer_stat_cor_plot(cor_stats_layer)
+#' 
+#' ## add colors
+#' ## add libd_layer_colors to refrence Human Pilot layers
+#' layer_stat_cor_plot(cor_stats_layer, reference_colors = libd_layer_colors)
+#' 
+#' ## supply polychrome colors to query clusters
+#' cluster_colors <-  c('#5A5156', '#E4E1E3', '#F6222E', '#FE00FA', '#16FF32', '#3283FE', '#FEAF16', '#B00068', '#1CFFCE')
+#' names(cluster_colors) <- rownames(cor_stats_layer)
 #'
-#' ## Repeat with just the top 10 layer marker genes
-#' layer_stat_cor_plot(layer_stat_cor(
-#'     tstats_Human_DLPFC_snRNAseq_Nguyen_topLayer,
-#'     modeling_results,
-#'     model_type = "enrichment",
-#'     top_n = 10
-#' ), max = 0.25)
-#'
-#' ## Now with the "pairwise" modeling results and also top_n = 10
-#' layer_stat_cor_plot(layer_stat_cor(
-#'     tstats_Human_DLPFC_snRNAseq_Nguyen_topLayer,
-#'     modeling_results,
-#'     model_type = "pairwise",
-#'     top_n = 10
-#' ), max = 0.25)
-layer_stat_cor_plot <-
-    function(
-        cor_stats_layer,
-        max = 0.81,
-        min = -max,
-        layerHeights = NULL,
-        cex = 1.2) {
-        ## From https://github.com/LieberInstitute/HumanPilot/blob/master/Analysis/Layer_Guesses/dlpfc_snRNAseq_annotation.R
-        theSeq <- seq(min, max, by = 0.01)
-        my.col <- grDevices::colorRampPalette(RColorBrewer::brewer.pal(7, "PRGn"))(length(theSeq))
-
-        ## Subset values
-        cor_stats_layer[cor_stats_layer <= min] <- min
-        cor_stats_layer[cor_stats_layer >= max] <- max
-
-        ## Re-shape the matrix
-        mat_vals <- t(cor_stats_layer)
-
-        ## Re-order and shorten names if they match our data
-        if (all(rownames(mat_vals) %in% c("WM", paste0("Layer", seq_len(6))))) {
-            rownames(mat_vals) <- gsub("ayer", "", rownames(mat_vals))
-            mat_vals <- mat_vals[c("WM", paste0("L", rev(seq_len(6)))), , drop = FALSE]
+#' layer_stat_cor_plot(cor_stats_layer, 
+#'                             query_colors = cluster_colors,
+#'                             reference_colors = libd_layer_colors)
+#' 
+#' ## Apply additional ComplexHeatmap param
+#' layer_stat_cor_plot(cor_stats_layer, cluster_rows = FALSE, cluster_columns = FALSE)
+#' 
+#' ## Add annotation
+#' annotation_df <- annotate_registered_clusters(cor_stats_layer, confidence_threshold = .55)
+#' layer_stat_cor_plot(cor_stats_layer, annotation = annotation_df)
+#' 
+#' ## All together
+#' layer_stat_cor_plot(cor_stats_layer, 
+#'                             query_colors = cluster_colors,
+#'                             reference_colors = libd_layer_colors, 
+#'                             annotation = annotation_df, 
+#'                             cluster_rows = FALSE, 
+#'                             cluster_columns = FALSE)
+#'                             
+layer_stat_cor_plot <- function(cor_stats_layer,
+                                        color_max = max(cor_stats_layer),
+                                        color_min = min(cor_stats_layer),
+                                        color_scale = RColorBrewer::brewer.pal(7, "PRGn"),
+                                        query_colors = NULL,
+                                        reference_colors = NULL,
+                                        annotation = NULL,
+                                        ...
+){
+
+  ## define color pallet
+  theSeq = seq(color_min, color_max, by = 0.01)
+  my.col = grDevices::colorRampPalette(color_scale)(length(theSeq))
+
+  # ## query annotations on row
+  if(!is.null(query_colors)){
+
+    stopifnot(all(rownames(cor_stats_layer) %in% names(query_colors)))
+    query_colors <- query_colors[rownames(cor_stats_layer)]
+
+
+    query_row_annotation <- ComplexHeatmap::rowAnnotation(
+      " " = rownames(cor_stats_layer),
+      col = list(" " = query_colors),
+      show_legend = FALSE)
+
+  } else query_row_annotation <- NULL
+
+  ## reference annotation on bottom
+  if(!is.null(reference_colors)){
+    stopifnot(all(colnames(cor_stats_layer) %in% names(reference_colors)))
+    reference_colors <- reference_colors[colnames(cor_stats_layer)]
+
+    ref_col_annotation <-  ComplexHeatmap::columnAnnotation(
+      " " = colnames(cor_stats_layer),
+      col = list(" " = reference_colors),
+      show_legend = FALSE
+    )
+  } else ref_col_annotation <- NULL
+
+  ## add annotation
+  if(!is.null(annotation)){
+    anno_matrix <- create_annotation_matrix(annotation, cor_stats_layer)
+
+    ## plot heatmap
+    return(
+      ComplexHeatmap::Heatmap(
+        matrix = cor_stats_layer,
+        col = my.col,
+        name = "Cor",
+        bottom_annotation = ref_col_annotation,
+        right_annotation = query_row_annotation,
+        cell_fun = function(j, i, x, y, width, height, fill) {
+          grid.text(anno_matrix[i, j], x, y, gp = gpar(fontsize = 10))
+        },
+        ...
+      ))
+  }
+
+  ## plot heatmap
+  return(
+    ComplexHeatmap::Heatmap(
+      matrix = cor_stats_layer,
+      col = my.col,
+      name = "Cor",
+      bottom_annotation = ref_col_annotation,
+      right_annotation = query_row_annotation,
+      ...
+    ))
+
+
+}
 
-            ## Use our default layer heights also
-            if (is.null(layerHeights)) {
-                layerHeights <- c(0, 40, 55, 75, 85, 110, 120, 135)
-            }
-        }
+create_annotation_matrix <- function(annotation_df, cor_stats_layer){
+
+  anno_list <- lapply(rownames(cor_stats_layer), 
+                      function(cluster){
+                        # look up confidence                    
+                        confidence <- annotation_df[match(cluster, annotation_df$cluster),"layer_confidence"]
+                        sym <- ifelse(confidence=="good", "X","*")
+                        # match annotations
+                        anno <- annotation_df[match(cluster, annotation_df$cluster),"layer_label"]
+                        return(ifelse(unlist(lapply(colnames(cor_stats_layer), grepl, anno)),sym,""))
+                      })
+
+  anno_matrix <- t(data.frame(anno_list))
+  rownames(anno_matrix) <- rownames(cor_stats_layer)
+  colnames(anno_matrix) <- colnames(cor_stats_layer)
+
+  return(anno_matrix)
+}
 
-        ## From fields:::imagePlotInfo
-        midpoints <- seq(min, max, length.out = length(my.col))
-        delta <- (midpoints[2] - midpoints[1]) / 2
-        breaks <- c(midpoints[1] - delta, midpoints + delta)
 
-        legend_cuts <- seq(-1, 1, by = 0.1)
-        legend_cuts <- legend_cuts[legend_cuts >= min & legend_cuts <= max]
-        axis.args <- list(
-            at = legend_cuts,
-            labels = legend_cuts
-        )
 
-        layer_matrix_plot(
-            matrix_values = mat_vals,
-            matrix_labels = NULL,
-            xlabs = NULL,
-            layerHeights = layerHeights,
-            mypal = my.col,
-            breaks = breaks,
-            axis.args = axis.args,
-            srt = 90,
-            cex = cex
-        )
-    }