diff --git a/DESCRIPTION b/DESCRIPTION
index a95b2b05..ad66f27f 100644
--- a/DESCRIPTION
+++ b/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
diff --git a/NAMESPACE b/NAMESPACE
index 42d8758a..c85f884a 100644
--- a/NAMESPACE
+++ b/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)
diff --git a/R/annotate_registered_clusters.R b/R/annotate_registered_clusters.R
index 13e1d5ab..865343c4 100644
--- a/R/annotate_registered_clusters.R
+++ b/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)
     }
 
diff --git a/R/layer_stat_cor_plot.R b/R/layer_stat_cor_plot.R
index 0d2653b8..df432d8a 100644
--- a/R/layer_stat_cor_plot.R
+++ b/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
-        )
-    }
diff --git a/man/layer_stat_cor_plot.Rd b/man/layer_stat_cor_plot.Rd
index 73496b45..8eccd53d 100644
--- a/man/layer_stat_cor_plot.Rd
+++ b/man/layer_stat_cor_plot.Rd
@@ -2,98 +2,115 @@
 % Please edit documentation in R/layer_stat_cor_plot.R
 \name{layer_stat_cor_plot}
 \alias{layer_stat_cor_plot}
-\title{Visualize the layer modeling correlation of statistics}
+\title{Visualize the correlation of layer modeling t-statistics with ComplexHeatmap}
 \usage{
 layer_stat_cor_plot(
   cor_stats_layer,
-  max = 0.81,
-  min = -max,
-  layerHeights = NULL,
-  cex = 1.2
+  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,
+  ...
 )
 }
 \arguments{
 \item{cor_stats_layer}{The output of \code{\link[=layer_stat_cor]{layer_stat_cor()}}.}
 
-\item{max}{A \code{numeric(1)} specifying the highest correlation value for the
-color scale (should be between 0 and 1).}
+\item{color_max}{A \code{numeric(1)} specifying the highest correlation value for
+the color scale (should be between 0 and 1).}
 
-\item{min}{A \code{numeric(1)} specifying the lowest correlation value for the
-color scale (should be between 0 and -1).}
+\item{color_min}{A \code{numeric(1)} specifying the lowest correlation value for
+the color scale (should be between 0 and -1).}
 
-\item{layerHeights}{A \code{numeric()} vector of length equal to
-\code{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 \code{\link[=layer_matrix_plot]{layer_matrix_plot()}}.}
+\item{color_scale}{A \code{character} vector specifying the color scale for the
+fill of the heatmap, defaults to classic purple -> green.}
 
-\item{cex}{Passed to \code{\link[=layer_matrix_plot]{layer_matrix_plot()}}.}
+\item{query_colors}{named \code{character} vector of colors, Adds colors to query
+row annotations.}
+
+\item{reference_colors}{named \code{character} vector of colors, Adds colors to
+reference column annotations.}
+
+\item{annotation}{annotation data.frame output of \code{\link[=annotate_registered_clusters]{annotate_registered_clusters()}},
+adds 'X' for good confidence annotations, '*' for poor confidence.}
+
+\item{...}{Additional parameters passed to \code{\link[ComplexHeatmap:Heatmap]{ComplexHeatmap::Heatmap()}}
+ex. \code{cluster_rows} and \code{cluster_columns}.}
 }
 \value{
-A heatmap for the correlation matrix between statistics.
+(\link[ComplexHeatmap:Heatmap-class]{Heatmap-class}) plot of t-stat correlations
 }
 \description{
-This function makes a heatmap from the \code{\link[=layer_stat_cor]{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).
+Use ComplexHeatmap to plot the correlation matrix between a reference and
+query modeling statistics from \code{\link[=layer_stat_cor]{layer_stat_cor()}}.
 }
 \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.
+Includes functionality to add color annotations,
+(helpful to match to colors in Visium spot plots), and annotations from
+\code{\link[=annotate_registered_clusters]{annotate_registered_clusters()}}.
 }
 \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))
+## most basic
+layer_stat_cor_plot(cor_stats_layer)
 
-## Restrict the range of colors further
-layer_stat_cor_plot(cor_stats_layer, max = 0.25)
+## add colors
+## add libd_layer_colors to refrence Human Pilot layers
+layer_stat_cor_plot(cor_stats_layer, reference_colors = libd_layer_colors)
 
-## 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)
+## 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)
 
-## 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(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)
+                            
 }
 \seealso{
-layer_matrix_plot annotate_registered_clusters
-
 Other Layer correlation functions: 
 \code{\link{annotate_registered_clusters}()},
 \code{\link{layer_stat_cor}()}
 }
 \author{
-Andrew E Jaffe, Leonardo Collado-Torres
+Louise Huuki-Myers
 }
 \concept{Layer correlation functions}