Multiple group bars on heatmap #2201

zehualilab · 2019-10-14T16:54:16Z

It seems that the Doheatmap function provides only one group bar coloring one metadata. Is there any way to add more group bars to color different metadata in one plot?

satijalab · 2019-10-27T16:40:21Z

We don't support this in the package, primarily because we would have to switch to other heatmap plotting functions which can be dramatically slower. For example, I believe the aheatmap (http://nmf.r-forge.r-project.org/aheatmap.html) package supports this, but I'm not quite sure. If you do find other packages that support this, please feel free to post here.

arkal · 2019-11-14T22:59:10Z

FWIW, I modified DoHeatmap to do this in a simple way. I didn't get to adding a legend for the additional groups but the bar on the top can be seen to better understand your data.

E.g.: MySeuratObj integrated 3 samples together and i wanted to see the DGE per identified cluster (seurat_idents) and additionally wanted to see the sample each column in the heatmap came from (orig.ident) so I can do

> DoMultiBarHeatmap(MySeuratObj, assay = 'SCT', features = top_10$gene, group.by='seurat_idents', additional.group.by = 'orig.ident')

This is the result (genes and clusters, etc are redacted since this is unpublished work ;) )

Here's the code

suppressPackageStartupMessages({
  library(rlang)
})

DoMultiBarHeatmap <- function (object, 
                               features = NULL, 
                               cells = NULL, 
                               group.by = "ident", 
                               additional.group.by = NULL, 
                               group.bar = TRUE, 
                               disp.min = -2.5, 
                               disp.max = NULL, 
                               slot = "scale.data", 
                               assay = NULL, 
                               label = TRUE, 
                               size = 5.5, 
                               hjust = 0, 
                               angle = 45, 
                               raster = TRUE, 
                               draw.lines = TRUE, 
                               lines.width = NULL, 
                               group.bar.height = 0.02, 
                               combine = TRUE) 
{
  cells <- cells %||% colnames(x = object)
  if (is.numeric(x = cells)) {
    cells <- colnames(x = object)[cells]
  }
  assay <- assay %||% DefaultAssay(object = object)
  DefaultAssay(object = object) <- assay
  features <- features %||% VariableFeatures(object = object)
  ## Why reverse???
  features <- rev(x = unique(x = features))
  disp.max <- disp.max %||% ifelse(test = slot == "scale.data", 
                                   yes = 2.5, no = 6)
  possible.features <- rownames(x = GetAssayData(object = object, 
                                                 slot = slot))
  if (any(!features %in% possible.features)) {
    bad.features <- features[!features %in% possible.features]
    features <- features[features %in% possible.features]
    if (length(x = features) == 0) {
      stop("No requested features found in the ", slot, 
           " slot for the ", assay, " assay.")
    }
    warning("The following features were omitted as they were not found in the ", 
            slot, " slot for the ", assay, " assay: ", paste(bad.features, 
                                                             collapse = ", "))
  }
  data <- as.data.frame(x = as.matrix(x = t(x = GetAssayData(object = object, 
                                                             slot = slot)[features, cells, drop = FALSE])))

  object <- suppressMessages(expr = StashIdent(object = object, 
                                               save.name = "ident"))
  group.by <- group.by %||% "ident"
  groups.use <- object[[c(group.by, additional.group.by)]][cells, , drop = FALSE]
  plots <- list()
  for (i in group.by) {
    data.group <- data
    group.use <- groups.use[, c(i, additional.group.by), drop = FALSE]
    
    for(colname in colnames(group.use)){
      if (!is.factor(x = group.use[[colname]])) {
        group.use[[colname]] <- factor(x = group.use[[colname]])
      }  
    }
    
    if (draw.lines) {
      lines.width <- lines.width %||% ceiling(x = nrow(x = data.group) * 
                                                0.0025)
      placeholder.cells <- sapply(X = 1:(length(x = levels(x = group.use[[i]])) * 
                                           lines.width), FUN = function(x) {
                                             return(Seurat:::RandomName(length = 20))
                                           })
      placeholder.groups <- data.frame(foo=rep(x = levels(x = group.use[[i]]), times = lines.width))
      placeholder.groups[additional.group.by] = NA
      colnames(placeholder.groups) <- colnames(group.use)
      rownames(placeholder.groups) <- placeholder.cells
      
      group.levels <- levels(x = group.use[[i]])

      group.use <- sapply(group.use, as.vector)
      rownames(x = group.use) <- cells
      
      group.use <- rbind(group.use, placeholder.groups)
      
      na.data.group <- matrix(data = NA, nrow = length(x = placeholder.cells), 
                              ncol = ncol(x = data.group), dimnames = list(placeholder.cells, 
                                                                           colnames(x = data.group)))
      data.group <- rbind(data.group, na.data.group)
    }
    
    #group.use = group.use[order(group.use[[i]]), , drop=F]
    group.use <- group.use[with(group.use, eval(parse(text=paste('order(', paste(c(i, additional.group.by), collapse=', '), ')', sep='')))), , drop=F]
    
    plot <- Seurat:::SingleRasterMap(data = data.group, raster = raster, 
                            disp.min = disp.min, disp.max = disp.max, feature.order = features, 
                            cell.order = rownames(x = group.use), group.by = group.use[[i]])
    
    if (group.bar) {
      pbuild <- ggplot_build(plot = plot)
      group.use2 <- group.use
      cols <- list()
      na.group <- Seurat:::RandomName(length = 20)
      for (colname in rev(x = colnames(group.use2))){
        if (colname == group.by){
          colid = paste0('Identity (', colname, ')')
        } else {
          colid = colname
        }
        
        if (draw.lines) {
          levels(x = group.use2[[colname]]) <- c(levels(x = group.use2[[colname]]), na.group)  
          group.use2[placeholder.cells, colname] <- na.group
          cols[[colname]] <- c(scales::hue_pal()(length(x = levels(x = group.use[[colname]]))), "#FFFFFF")
        } else {
          cols[[colname]] <- c(scales::hue_pal()(length(x = levels(x = group.use[[colname]]))))
        }
        names(x = cols[[colname]]) <- levels(x = group.use2[[colname]])
        
        
        y.range <- diff(x = pbuild$layout$panel_params[[1]]$y.range)
        y.pos <- max(pbuild$layout$panel_params[[1]]$y.range) + y.range * 0.015
        y.max <- y.pos + group.bar.height * y.range
        pbuild$layout$panel_params[[1]]$y.range <- c(pbuild$layout$panel_params[[1]]$y.range[1], y.max)
        
        plot <- suppressMessages(plot + 
                                annotation_raster(raster = t(x = cols[[colname]][group.use2[[colname]]]),  xmin = -Inf, xmax = Inf, ymin = y.pos, ymax = y.max) + 
                                annotation_custom(grob = grid::textGrob(label = colid, hjust = 0, gp = grid::gpar(cex = 0.75)), ymin = mean(c(y.pos, y.max)), ymax = mean(c(y.pos, y.max)), xmin = Inf, xmax = Inf) +
                                coord_cartesian(ylim = c(0, y.max), clip = "off")) 
        
        #temp <- as.data.frame(cols[[colname]][levels(group.use[[colname]])])
        #colnames(temp) <- 'color'
        #temp$x <- temp$y <- 1
        #temp[['name']] <- as.factor(rownames(temp))
        
        #temp <- ggplot(temp, aes(x=x, y=y, fill=name)) + geom_point(shape=21, size=5) + labs(fill=colname) + theme(legend.position = "bottom")
        #legend <- get_legend(temp)
        #multiplot(plot, legend, heights=3,1)
        
        if ((colname == group.by) && label) {
          x.max <- max(pbuild$layout$panel_params[[1]]$x.range)
          x.divs <- pbuild$layout$panel_params[[1]]$x.major
          group.use$x <- x.divs
          label.x.pos <- tapply(X = group.use$x, INDEX = group.use[[colname]],
                                FUN = median) * x.max
          label.x.pos <- data.frame(group = names(x = label.x.pos), 
                                    label.x.pos)
          plot <- plot + geom_text(stat = "identity", 
                                   data = label.x.pos, aes_string(label = "group", 
                                                                  x = "label.x.pos"), y = y.max + y.max * 
                                     0.03 * 0.5, angle = angle, hjust = hjust, 
                                   size = size)
          plot <- suppressMessages(plot + coord_cartesian(ylim = c(0, 
                                                                   y.max + y.max * 0.002 * max(nchar(x = levels(x = group.use[[colname]]))) * 
                                                                     size), clip = "off"))
        }
      }
    }
    plot <- plot + theme(line = element_blank())
    plots[[i]] <- plot
  }
  if (combine) {
    plots <- CombinePlots(plots = plots)
  }
  return(plots)
}

HomairaH · 2019-12-13T17:51:54Z

@arkal Thanks for this code, very helpful! I have used it on my data to add two colour bars on the heatmap, however I am having trouble with two aspects:

Let's say I have clusters 1, 2, 3, 4, 5. I reordered them to 2, 1, 3, 5, 4 and added this order as a column in the metadata called "clusters_new".

However, when I change the Idents of the object to "clusters_new" and use "clusters_new" in the group.by section, the order of the heatmap still appears as 1, 2, 3, 4, 5 instead of 2, 1, 3, 5, 4. Is there a way to switch the order the clusters appear on the heatmap?

Is it possible to specify the colors for the color bar for both the group.by and additional.group.by clusters?

FWIW, I modified DoHeatmap to do this in a simple way. I didn't get to adding a legend for the additional groups but the bar on the top can be seen to better understand your data.

E.g.: MySeuratObj integrated 3 samples together and i wanted to see the DGE per identified cluster (seurat_idents) and additionally wanted to see the sample each column in the heatmap came from (orig.ident) so I can do

> DoMultiBarHeatmap(MySeuratObj, assay = 'SCT', features = top_10$gene, group.by='seurat_idents', additional.group.by = 'orig.ident')

This is the result (genes and clusters, etc are redacted since this is unpublished work ;) )

Here's the code

suppressPackageStartupMessages({
  library(rlang)
})

DoMultiBarHeatmap <- function (object, 
                               features = NULL, 
                               cells = NULL, 
                               group.by = "ident", 
                               additional.group.by = NULL, 
                               group.bar = TRUE, 
                               disp.min = -2.5, 
                               disp.max = NULL, 
                               slot = "scale.data", 
                               assay = NULL, 
                               label = TRUE, 
                               size = 5.5, 
                               hjust = 0, 
                               angle = 45, 
                               raster = TRUE, 
                               draw.lines = TRUE, 
                               lines.width = NULL, 
                               group.bar.height = 0.02, 
                               combine = TRUE) 
{
  cells <- cells %||% colnames(x = object)
  if (is.numeric(x = cells)) {
    cells <- colnames(x = object)[cells]
  }
  assay <- assay %||% DefaultAssay(object = object)
  DefaultAssay(object = object) <- assay
  features <- features %||% VariableFeatures(object = object)
  ## Why reverse???
  features <- rev(x = unique(x = features))
  disp.max <- disp.max %||% ifelse(test = slot == "scale.data", 
                                   yes = 2.5, no = 6)
  possible.features <- rownames(x = GetAssayData(object = object, 
                                                 slot = slot))
  if (any(!features %in% possible.features)) {
    bad.features <- features[!features %in% possible.features]
    features <- features[features %in% possible.features]
    if (length(x = features) == 0) {
      stop("No requested features found in the ", slot, 
           " slot for the ", assay, " assay.")
    }
    warning("The following features were omitted as they were not found in the ", 
            slot, " slot for the ", assay, " assay: ", paste(bad.features, 
                                                             collapse = ", "))
  }
  data <- as.data.frame(x = as.matrix(x = t(x = GetAssayData(object = object, 
                                                             slot = slot)[features, cells, drop = FALSE])))

  object <- suppressMessages(expr = StashIdent(object = object, 
                                               save.name = "ident"))
  group.by <- group.by %||% "ident"
  groups.use <- object[[c(group.by, additional.group.by)]][cells, , drop = FALSE]
  plots <- list()
  for (i in group.by) {
    data.group <- data
    group.use <- groups.use[, c(i, additional.group.by), drop = FALSE]
    
    for(colname in colnames(group.use)){
      if (!is.factor(x = group.use[[colname]])) {
        group.use[[colname]] <- factor(x = group.use[[colname]])
      }  
    }
    
    if (draw.lines) {
      lines.width <- lines.width %||% ceiling(x = nrow(x = data.group) * 
                                                0.0025)
      placeholder.cells <- sapply(X = 1:(length(x = levels(x = group.use[[i]])) * 
                                           lines.width), FUN = function(x) {
                                             return(Seurat:::RandomName(length = 20))
                                           })
      placeholder.groups <- data.frame(foo=rep(x = levels(x = group.use[[i]]), times = lines.width))
      placeholder.groups[additional.group.by] = NA
      colnames(placeholder.groups) <- colnames(group.use)
      rownames(placeholder.groups) <- placeholder.cells
      
      group.levels <- levels(x = group.use[[i]])

      group.use <- sapply(group.use, as.vector)
      rownames(x = group.use) <- cells
      
      group.use <- rbind(group.use, placeholder.groups)
      
      na.data.group <- matrix(data = NA, nrow = length(x = placeholder.cells), 
                              ncol = ncol(x = data.group), dimnames = list(placeholder.cells, 
                                                                           colnames(x = data.group)))
      data.group <- rbind(data.group, na.data.group)
    }
    
    #group.use = group.use[order(group.use[[i]]), , drop=F]
    group.use <- group.use[with(group.use, eval(parse(text=paste('order(', paste(c(i, additional.group.by), collapse=', '), ')', sep='')))), , drop=F]
    
    plot <- Seurat:::SingleRasterMap(data = data.group, raster = raster, 
                            disp.min = disp.min, disp.max = disp.max, feature.order = features, 
                            cell.order = rownames(x = group.use), group.by = group.use[[i]])
    
    if (group.bar) {
      pbuild <- ggplot_build(plot = plot)
      group.use2 <- group.use
      cols <- list()
      na.group <- Seurat:::RandomName(length = 20)
      for (colname in rev(x = colnames(group.use2))){
        if (colname == group.by){
          colid = paste0('Identity (', colname, ')')
        } else {
          colid = colname
        }
        
        if (draw.lines) {
          levels(x = group.use2[[colname]]) <- c(levels(x = group.use2[[colname]]), na.group)  
          group.use2[placeholder.cells, colname] <- na.group
          cols[[colname]] <- c(scales::hue_pal()(length(x = levels(x = group.use[[colname]]))), "#FFFFFF")
        } else {
          cols[[colname]] <- c(scales::hue_pal()(length(x = levels(x = group.use[[colname]]))))
        }
        names(x = cols[[colname]]) <- levels(x = group.use2[[colname]])
        
        
        y.range <- diff(x = pbuild$layout$panel_params[[1]]$y.range)
        y.pos <- max(pbuild$layout$panel_params[[1]]$y.range) + y.range * 0.015
        y.max <- y.pos + group.bar.height * y.range
        pbuild$layout$panel_params[[1]]$y.range <- c(pbuild$layout$panel_params[[1]]$y.range[1], y.max)
        
        plot <- suppressMessages(plot + 
                                annotation_raster(raster = t(x = cols[[colname]][group.use2[[colname]]]),  xmin = -Inf, xmax = Inf, ymin = y.pos, ymax = y.max) + 
                                annotation_custom(grob = grid::textGrob(label = colid, hjust = 0, gp = grid::gpar(cex = 0.75)), ymin = mean(c(y.pos, y.max)), ymax = mean(c(y.pos, y.max)), xmin = Inf, xmax = Inf) +
                                coord_cartesian(ylim = c(0, y.max), clip = "off")) 
        
        #temp <- as.data.frame(cols[[colname]][levels(group.use[[colname]])])
        #colnames(temp) <- 'color'
        #temp$x <- temp$y <- 1
        #temp[['name']] <- as.factor(rownames(temp))
        
        #temp <- ggplot(temp, aes(x=x, y=y, fill=name)) + geom_point(shape=21, size=5) + labs(fill=colname) + theme(legend.position = "bottom")
        #legend <- get_legend(temp)
        #multiplot(plot, legend, heights=3,1)
        
        if ((colname == group.by) && label) {
          x.max <- max(pbuild$layout$panel_params[[1]]$x.range)
          x.divs <- pbuild$layout$panel_params[[1]]$x.major
          group.use$x <- x.divs
          label.x.pos <- tapply(X = group.use$x, INDEX = group.use[[colname]],
                                FUN = median) * x.max
          label.x.pos <- data.frame(group = names(x = label.x.pos), 
                                    label.x.pos)
          plot <- plot + geom_text(stat = "identity", 
                                   data = label.x.pos, aes_string(label = "group", 
                                                                  x = "label.x.pos"), y = y.max + y.max * 
                                     0.03 * 0.5, angle = angle, hjust = hjust, 
                                   size = size)
          plot <- suppressMessages(plot + coord_cartesian(ylim = c(0, 
                                                                   y.max + y.max * 0.002 * max(nchar(x = levels(x = group.use[[colname]]))) * 
                                                                     size), clip = "off"))
        }
      }
    }
    plot <- plot + theme(line = element_blank())
    plots[[i]] <- plot
  }
  if (combine) {
    plots <- CombinePlots(plots = plots)
  }
  return(plots)
}

arkal · 2019-12-17T01:55:37Z

@HomairaH I'm glad it helped you. Silly me I was recalculating levels instead of inheriting. I modified the code and The Code is at the bottom.

I added a new parameter additional.group.sort.by That allows you to specify that you'd like to sort cells additionally by groups in the new bar annotation.
I also added a parameter cols.use that takes a list of (optionally named) vectors of colors so you can provide colors to your plot.

#Order is G1, G2M, S
DoMultiBarHeatmap(sobj, features=top2$gene, group.by='Phase', additional.group.by = c('seurat_clusters', 'foo', 'bar'))

# Change order to S, G1, G2M
sobj@meta.data$Phase <- factor(x = sobj@meta.data$Phase, levels=c('S', 'G1', 'G2M'))
DoMultiBarHeatmap(sobj, features=top2$gene, group.by='Phase', additional.group.by = c('seurat_clusters', 'foo', 'bar'))

# Show we can sort sub-bars
DoMultiBarHeatmap(sobj, features=top2$gene, group.by='Phase', additional.group.by = c('seurat_clusters', 'foo', 'bar'), additional.group.sort.by = c('seurat_clusters'))

# Show we can color anything
cols.use <- list(Phase=c('red', 'blue', 'pink', 'brown'))
DoMultiBarHeatmap(sobj, features=top2$gene, group.by='Phase', additional.group.by = c('seurat_clusters', 'foo', 'bar'), additional.group.sort.by = c('bar'), cols.use=cols.use)

names(cols.use[['Phase']]) <- c('S', 'foo', 'G2M', 'G1')
DoMultiBarHeatmap(sobj, features=top2$gene, group.by='Phase', additional.group.by = c('seurat_clusters', 'foo', 'bar'), additional.group.sort.by = c('bar'), cols.use=cols.use)

Here's all the code!

suppressPackageStartupMessages({
  library(rlang)
})

DoMultiBarHeatmap <- function (object, 
                               features = NULL, 
                               cells = NULL, 
                               group.by = "ident", 
                               additional.group.by = NULL, 
                               additional.group.sort.by = NULL, 
                               cols.use = NULL,
                               group.bar = TRUE, 
                               disp.min = -2.5, 
                               disp.max = NULL, 
                               slot = "scale.data", 
                               assay = NULL, 
                               label = TRUE, 
                               size = 5.5, 
                               hjust = 0, 
                               angle = 45, 
                               raster = TRUE, 
                               draw.lines = TRUE, 
                               lines.width = NULL, 
                               group.bar.height = 0.02, 
                               combine = TRUE) 
{
  cells <- cells %||% colnames(x = object)
  if (is.numeric(x = cells)) {
    cells <- colnames(x = object)[cells]
  }
  assay <- assay %||% DefaultAssay(object = object)
  DefaultAssay(object = object) <- assay
  features <- features %||% VariableFeatures(object = object)
  ## Why reverse???
  features <- rev(x = unique(x = features))
  disp.max <- disp.max %||% ifelse(test = slot == "scale.data", 
                                   yes = 2.5, no = 6)
  possible.features <- rownames(x = GetAssayData(object = object, 
                                                 slot = slot))
  if (any(!features %in% possible.features)) {
    bad.features <- features[!features %in% possible.features]
    features <- features[features %in% possible.features]
    if (length(x = features) == 0) {
      stop("No requested features found in the ", slot, 
           " slot for the ", assay, " assay.")
    }
    warning("The following features were omitted as they were not found in the ", 
            slot, " slot for the ", assay, " assay: ", paste(bad.features, 
                                                             collapse = ", "))
  }
  
  if (!is.null(additional.group.sort.by)) {
    if (any(!additional.group.sort.by %in% additional.group.by)) {
      bad.sorts <- additional.group.sort.by[!additional.group.sort.by %in% additional.group.by]
      additional.group.sort.by <- additional.group.sort.by[additional.group.sort.by %in% additional.group.by]
      if (length(x = bad.sorts) > 0) {
        warning("The following additional sorts were omitted as they were not a subset of additional.group.by : ", 
                paste(bad.sorts, collapse = ", "))
      }
    }
  }
  
  data <- as.data.frame(x = as.matrix(x = t(x = GetAssayData(object = object, 
                                                             slot = slot)[features, cells, drop = FALSE])))
  
  object <- suppressMessages(expr = StashIdent(object = object, 
                                               save.name = "ident"))
  group.by <- group.by %||% "ident"
  groups.use <- object[[c(group.by, additional.group.by[!additional.group.by %in% group.by])]][cells, , drop = FALSE]
  plots <- list()
  for (i in group.by) {
    data.group <- data
    if (!is_null(additional.group.by)) {
      additional.group.use <- additional.group.by[additional.group.by!=i]  
      if (!is_null(additional.group.sort.by)){
        additional.sort.use = additional.group.sort.by[additional.group.sort.by != i]  
      } else {
        additional.sort.use = NULL
      }
    } else {
      additional.group.use = NULL
      additional.sort.use = NULL
    }
    
    group.use <- groups.use[, c(i, additional.group.use), drop = FALSE]
    
    for(colname in colnames(group.use)){
      if (!is.factor(x = group.use[[colname]])) {
        group.use[[colname]] <- factor(x = group.use[[colname]])
      }  
    }
    
    if (draw.lines) {
      lines.width <- lines.width %||% ceiling(x = nrow(x = data.group) * 
                                                0.0025)
      placeholder.cells <- sapply(X = 1:(length(x = levels(x = group.use[[i]])) * 
                                           lines.width), FUN = function(x) {
                                             return(Seurat:::RandomName(length = 20))
                                           })
      placeholder.groups <- data.frame(rep(x = levels(x = group.use[[i]]), times = lines.width))
      group.levels <- list()
      group.levels[[i]] = levels(x = group.use[[i]])
      for (j in additional.group.use) {
        group.levels[[j]] <- levels(x = group.use[[j]])
        placeholder.groups[[j]] = NA
      }
      
      colnames(placeholder.groups) <- colnames(group.use)
      rownames(placeholder.groups) <- placeholder.cells
      
      group.use <- sapply(group.use, as.vector)
      rownames(x = group.use) <- cells
      
      group.use <- rbind(group.use, placeholder.groups)
      
      for (j in names(group.levels)) {
        group.use[[j]] <- factor(x = group.use[[j]], levels = group.levels[[j]])
      }
      
      na.data.group <- matrix(data = NA, nrow = length(x = placeholder.cells), 
                              ncol = ncol(x = data.group), dimnames = list(placeholder.cells, 
                                                                           colnames(x = data.group)))
      data.group <- rbind(data.group, na.data.group)
    }
    
    order_expr <- paste0('order(', paste(c(i, additional.sort.use), collapse=','), ')')
    group.use = with(group.use, group.use[eval(parse(text=order_expr)), , drop=F])
    
    plot <- Seurat:::SingleRasterMap(data = data.group, raster = raster, 
                                     disp.min = disp.min, disp.max = disp.max, feature.order = features, 
                                     cell.order = rownames(x = group.use), group.by = group.use[[i]])
    
    if (group.bar) {
      pbuild <- ggplot_build(plot = plot)
      group.use2 <- group.use
      cols <- list()
      na.group <- Seurat:::RandomName(length = 20)
      for (colname in rev(x = colnames(group.use2))) {
        if (colname == i) {
          colid = paste0('Identity (', colname, ')')
        } else {
          colid = colname
        }
        
        # Default
        cols[[colname]] <- c(scales::hue_pal()(length(x = levels(x = group.use[[colname]]))))  
        
        #Overwrite if better value is provided
        if (!is_null(cols.use[[colname]])) {
          req_length = length(x = levels(group.use))
          if (length(cols.use[[colname]]) < req_length){
            warning("Cannot use provided colors for ", colname, " since there aren't enough colors.")
          } else {
            if (!is_null(names(cols.use[[colname]]))) {
              if (all(levels(group.use[[colname]]) %in% names(cols.use[[colname]]))) {
                cols[[colname]] <- as.vector(cols.use[[colname]][levels(group.use[[colname]])])
              } else {
                warning("Cannot use provided colors for ", colname, " since all levels (", paste(levels(group.use[[colname]]), collapse=","), ") are not represented.")
              }
            } else {
              cols[[colname]] <- as.vector(cols.use[[colname]])[c(1:length(x = levels(x = group.use[[colname]])))]
            }
          }
        }
        
        # Add white if there's lines
        if (draw.lines) {
          levels(x = group.use2[[colname]]) <- c(levels(x = group.use2[[colname]]), na.group)  
          group.use2[placeholder.cells, colname] <- na.group
          cols[[colname]] <- c(cols[[colname]], "#FFFFFF")
        }
        names(x = cols[[colname]]) <- levels(x = group.use2[[colname]])
        
        y.range <- diff(x = pbuild$layout$panel_params[[1]]$y.range)
        y.pos <- max(pbuild$layout$panel_params[[1]]$y.range) + y.range * 0.015
        y.max <- y.pos + group.bar.height * y.range
        pbuild$layout$panel_params[[1]]$y.range <- c(pbuild$layout$panel_params[[1]]$y.range[1], y.max)
        
        plot <- suppressMessages(plot + 
                                   annotation_raster(raster = t(x = cols[[colname]][group.use2[[colname]]]),  xmin = -Inf, xmax = Inf, ymin = y.pos, ymax = y.max) + 
                                   annotation_custom(grob = grid::textGrob(label = colid, hjust = 0, gp = gpar(cex = 0.75)), ymin = mean(c(y.pos, y.max)), ymax = mean(c(y.pos, y.max)), xmin = Inf, xmax = Inf) +
                                   coord_cartesian(ylim = c(0, y.max), clip = "off")) 
        
        if ((colname == i) && label) {
          x.max <- max(pbuild$layout$panel_params[[1]]$x.range)
          x.divs <- pbuild$layout$panel_params[[1]]$x.major
          group.use$x <- x.divs
          label.x.pos <- tapply(X = group.use$x, INDEX = group.use[[colname]],
                                FUN = median) * x.max
          label.x.pos <- data.frame(group = names(x = label.x.pos), 
                                    label.x.pos)
          plot <- plot + geom_text(stat = "identity", 
                                   data = label.x.pos, aes_string(label = "group", 
                                                                  x = "label.x.pos"), y = y.max + y.max * 
                                     0.03 * 0.5, angle = angle, hjust = hjust, 
                                   size = size)
          plot <- suppressMessages(plot + coord_cartesian(ylim = c(0, 
                                                                   y.max + y.max * 0.002 * max(nchar(x = levels(x = group.use[[colname]]))) * 
                                                                     size), clip = "off"))
        }
      }
    }
    plot <- plot + theme(line = element_blank())
    plots[[i]] <- plot
  }
  if (combine) {
    plots <- CombinePlots(plots = plots)
  }
  return(plots)
}

HomairaH · 2019-12-29T23:50:25Z

@arkal thanks so much again!

MinBio · 2020-01-15T22:16:34Z

Thanks for the code @arkal . It's very helpful. I noticed that the legend is missing for the new group bars. Could you add it to the heatmap? thank you.

arkal · 2020-01-16T22:22:41Z

Yeah that's what i mentioned upfront, it's a feature i don't have the time to add. Feel free to modify that code to do it :) FWIW the new group bars are ordered the same as your metadata object.

Maybe i'll have free time in the future to get to this but not right now sorry :(

kaizen89 · 2020-03-26T21:33:33Z

Hi @arkal , I was using your function DoMultiBarHeatmap but recently after updating some packages, it is not working anymore. It throws this error

Error in tapply(X = group.use$x, INDEX = group.use[[colname]], FUN = median) : 
  arguments must have same length

Do you know what might be causing this error? Thanks

ktessema · 2020-03-27T09:36:30Z

Hey @kaizen89 - I was wondering the same! I'm a complete newbie, so not completely sure if this is the correct answer, but replacing one line of code did the trick for me. Looks like it's an issue with the ggplot2 version. I found the replacement line in the DoHeatmap code. You could also set label=FALSE when calling the function, if you don't care about that.

Replace: x.divs <- pbuild$layout$panel_params[[1]]$x.major

With: x.divs <- pbuild$layout$panel_params[[1]]$x.major %||% pbuild$layout$panel_params[[1]]$x$break_positions()

Hope this helps!

kaizen89 · 2020-03-27T12:45:31Z

Indeed I was setting label=F as a workaround, but now your solution fixes the issue, thanks @iheartfoosball.

csung331 · 2020-07-16T00:43:51Z

@arkal
I am very new to R coding. So apologies in advance if this sounds too basic.
"DoMultiBarHeatmap" function is exactly what I need for my data analysis. But I can't seem to make my R recognize the code. Below is the error message that I keep getting.
ㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡ
Error in DoMultiBarHeatmap(CtrlCis, assay = "RNA", features = top10$gene, :
could not find function "DoMultiBarHeatmap"
ㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡ
I thought I was missing a package. I installed "rlang", "pheatmap", "ggplot2", "cowplot", "dplyr", "patchwork", and "Seurat". And it still seemed like R wasn't recognizing this code. Is there a package that I need to install that I am missing.

Thanks!!

ktessema · 2020-07-16T00:46:02Z

@csung331 Did you run @arkal's code (end of 12/16/19 post) to load the function into your R session?

csung331 · 2020-07-16T00:49:15Z

@iheartfoosball Below is the code that I ran.
DoMultiBarHeatmap(SeuratObject, assay = 'RNA', features=top10$gene, group.by='clusters', additional.group.by = 'orig.identity')

ktessema · 2020-07-16T01:38:26Z

@csung331 Looks like you haven't loaded the function yet, so R does not recognize it when you try to call it. I'd suggest copying the code at the end of arkal's Dec 16 post (after the line that reads "Here's all the code") and pasting it into an empty R Script file. Then you can run that code to load the function into the R session. You can then save that file separately so that you can use it to load the function in the future (using the 'source' function, as described here).

csung331 · 2020-07-16T15:24:24Z

@iheartfoosball Great! Thank you. It seems like the code is working. But now I am confronted with a different issue which might be due to the way how I processed and integrated my datasets. Below is the error message that I'm getting.
ㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡ
Error in DoMultiBarHeatmap(SeuratOjbect, assay = "RNA", features = top10$gene, :
No requested features found in the scale.data slot for the RNA assay.
ㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡ
I believe with the current Seurat data analysis flow, SCTransformed data cannot be used , otherwise I get the above error message. I would have to use NormalizeData() to normalize the data for integration.

ktessema · 2020-07-16T17:22:20Z

@csung331 No prob!

If you want to use the values resulting from SCTransform, you can set the assay to 'SCT' instead of 'RNA' when you call the DoMultiBarHeatmap function. I believe you can also specify the slot when you call the function (eg assay='SCT', slot='data' or assay='SCT', slot='scale.data'). So you can decide which values you wish to plot and specify the relevant parameters when calling the function.

csung331 · 2020-07-16T18:56:01Z

@iheartfoosball Thanks for all your insight.

I SCT processed my data from individual group but used NormalizeData() when integrating all datasets. So I should use assay='RNA', slot='data'. The code that I would use is
ㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡ
DoMultiBarHeatmap(SeuratObject, assay='RNA', slot='data', features=top10$gene, group.by='orig.ident', additional.group.by = 'clusters')
ㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡ

But when I run this code, I also get the below error message.
ㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡ
Error in t.default(x = GetAssayData(object = object, slot = slot)[features, :
argument is not a matrix
ㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡ
Would you know how to correct this error?

ktessema · 2020-07-16T19:08:20Z

@csung331 no prob

hmm have you tried without the features argument? i would try with just the first 3 arguments (object, assay, slot) and see if it runs.

PrashINRA · 2020-07-29T15:40:58Z

Thanks for the code @arkal . It's very helpful. I noticed that the legend is missing for the new group bars. Could you add it to the heatmap? thank you.

Indeed, this is what I was looking for too, did you manage it?

Diozao666 · 2020-08-02T20:49:30Z

Thanks Arkal you are doing god's work

However, I got the warning message

` Error in gpar(cex = 0.75) : could not find function "gpar"

`
Is this from another package?

EDIT: nvm
I loaded the library grid and it worked, also did the iheartfoosball modification for it to work
thanks guys, saved me hours of trouble

PrashINRA · 2020-08-02T21:20:11Z

Thanks Arkal you are doing god's work

However, I got the warning message

` Error in gpar(cex = 0.75) : could not find function "gpar"

`
Is this from another package?

EDIT: nvm
I loaded the library grid and it worked, also did the iheartfoosball modification for it to work
thanks guys, saved me hours of trouble

I think you just have to load library(grid) to use the function gpar.

shaln · 2021-02-02T01:51:36Z

Has anyone tried running this recently? I keep getting the following error and have no idea how to fix it. Error still appears even without cols.use. I've included my sessionInfo() below in case it's helpful! Would appreciate any help I can get with this. Thanks!

cols.use <- list(condition=c("tomato", "dodgerblue"))
DoMultiBarHeatmap(subset(dataset, downsample=50), group.by = "condition", additional.group.by = "cell.identity", cols.use)

Error: Must request at least one colour from a hue palette.

> sessionInfo()
R version 4.0.3 (2020-10-10)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19041)

Matrix products: default

locale:
[1] LC_COLLATE=English_United Kingdom.1252  LC_CTYPE=English_United Kingdom.1252   
[3] LC_MONETARY=English_United Kingdom.1252 LC_NUMERIC=C                           
[5] LC_TIME=English_United Kingdom.1252    

attached base packages:
 [1] grid      parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] rlang_0.4.10                plyr_1.8.6                  dplyr_1.0.3                
 [4] ggpubr_0.4.0                SingleR_1.4.0               scran_1.18.3               
 [7] batchelor_1.6.2             scater_1.18.3               ggplot2_3.3.3              
[10] Seurat_3.2.3                SingleCellExperiment_1.12.0 SummarizedExperiment_1.20.0
[13] Biobase_2.50.0              GenomicRanges_1.42.0        GenomeInfoDb_1.26.2        
[16] IRanges_2.24.1              S4Vectors_0.28.1            BiocGenerics_0.36.0        
[19] MatrixGenerics_1.2.0        matrixStats_0.57.0          Matrix_1.3-2  
loaded via a namespace (and not attached):
  [1] readxl_1.3.1              backports_1.2.1           igraph_1.2.6             
  [4] lazyeval_0.2.2            splines_4.0.3             BiocParallel_1.24.1      
  [7] listenv_0.8.0             scattermore_0.7           digest_0.6.27            
 [10] htmltools_0.5.1.1         viridis_0.5.1             magrittr_2.0.1           
 [13] tensor_1.5                cluster_2.1.0             ROCR_1.0-11              
 [16] openxlsx_4.2.3            limma_3.46.0              globals_0.14.0           
 [19] colorspace_2.0-0          ggrepel_0.9.1             haven_2.3.1              
 [22] crayon_1.3.4              RCurl_1.98-1.2            jsonlite_1.7.2           
 [25] spatstat_1.64-1           spatstat.data_1.7-0       survival_3.2-7           
 [28] zoo_1.8-8                 glue_1.4.2                polyclip_1.10-0          
 [31] gtable_0.3.0              zlibbioc_1.36.0           XVector_0.30.0           
 [34] leiden_0.3.6              DelayedArray_0.16.1       car_3.0-10               
 [37] BiocSingular_1.6.0        future.apply_1.7.0        abind_1.4-5              
 [40] scales_1.1.1              edgeR_3.32.1              rstatix_0.6.0            
 [43] miniUI_0.1.1.1            Rcpp_1.0.6                viridisLite_0.3.0        
 [46] xtable_1.8-4              reticulate_1.18           dqrng_0.2.1              
 [49] foreign_0.8-81            rsvd_1.0.3                ResidualMatrix_1.0.0     
 [52] htmlwidgets_1.5.3         httr_1.4.2                RColorBrewer_1.1-2       
 [55] ellipsis_0.3.1            ica_1.0-2                 farver_2.0.3             
 [58] pkgconfig_2.0.3           scuttle_1.0.4             uwot_0.1.10              
 [61] deldir_0.2-9              locfit_1.5-9.4            labeling_0.4.2           
 [64] tidyselect_1.1.0          reshape2_1.4.4            later_1.1.0.1            
 [67] cellranger_1.1.0          munsell_0.5.0             tools_4.0.3              
 [70] generics_0.1.0            broom_0.7.3               ggridges_0.5.3           
 [73] stringr_1.4.0             fastmap_1.0.1             goftest_1.2-2            
 [76] fitdistrplus_1.1-3        zip_2.1.1                 purrr_0.3.4              
 [79] RANN_2.6.1                pbapply_1.4-3             future_1.21.0            
 [82] nlme_3.1-151              sparseMatrixStats_1.2.0   mime_0.9                 
 [85] compiler_4.0.3            rstudioapi_0.13           curl_4.3                 
 [88] beeswarm_0.2.3            plotly_4.9.3              png_0.1-7                
 [91] ggsignif_0.6.0            spatstat.utils_2.0-0      tibble_3.0.5             
 [94] statmod_1.4.35            stringi_1.5.3             forcats_0.5.0            
 [97] lattice_0.20-41           bluster_1.0.0             vctrs_0.3.6              
[100] pillar_1.4.7              lifecycle_0.2.0           lmtest_0.9-38            
[103] RcppAnnoy_0.0.18          BiocNeighbors_1.8.2       data.table_1.13.6        
[106] cowplot_1.1.1             bitops_1.0-6              irlba_2.3.3              
[109] httpuv_1.5.5              patchwork_1.1.1           R6_2.5.0                 
[112] promises_1.1.1            rio_0.5.16                KernSmooth_2.23-18       
[115] gridExtra_2.3             vipor_0.4.5               parallelly_1.23.0        
[118] codetools_0.2-18          MASS_7.3-53               withr_2.4.0              
[121] sctransform_0.3.2         GenomeInfoDbData_1.2.4    hms_1.0.0                
[124] mgcv_1.8-33               rpart_4.1-15              beachmat_2.6.4           
[127] tidyr_1.1.2               DelayedMatrixStats_1.12.2 carData_3.0-4            
[130] Rtsne_0.15                shiny_1.5.0               ggbeeswarm_0.6.0

ktessema · 2021-02-04T01:05:35Z

Has anyone tried running this recently? I keep getting the following error and have no idea how to fix it. Error still appears even without cols.use. I've included my sessionInfo() below in case it's helpful! Would appreciate any help I can get with this. Thanks!
cols.use <- list(condition=c("tomato", "dodgerblue"))
DoMultiBarHeatmap(subset(dataset, downsample=50), group.by = "condition", additional.group.by = "cell.identity", cols.use)

Error: Must request at least one colour from a hue palette.

Do you know if all your cells have both a 'condition' and 'cell.identity' value? Maybe check by doing:

unique(dataset$condition)
unique(dataset$cell.identity)

InMyGenes · 2021-02-11T20:38:17Z

Hello, whenever I run the code above from

FWIW, I modified DoHeatmap to do this in a simple way. I didn't get to adding a legend for the additional groups but the bar on the top can be seen to better understand your data.

E.g.: MySeuratObj integrated 3 samples together and i wanted to see the DGE per identified cluster (seurat_idents) and additionally wanted to see the sample each column in the heatmap came from (orig.ident) so I can do

> DoMultiBarHeatmap(MySeuratObj, assay = 'SCT', features = top_10$gene, group.by='seurat_idents', additional.group.by = 'orig.ident')

This is the result (genes and clusters, etc are redacted since this is unpublished work ;) )

Here's the code

suppressPackageStartupMessages({
  library(rlang)
})

DoMultiBarHeatmap <- function (object, 
                               features = NULL, 
                               cells = NULL, 
                               group.by = "ident", 
                               additional.group.by = NULL, 
                               group.bar = TRUE, 
                               disp.min = -2.5, 
                               disp.max = NULL, 
                               slot = "scale.data", 
                               assay = NULL, 
                               label = TRUE, 
                               size = 5.5, 
                               hjust = 0, 
                               angle = 45, 
                               raster = TRUE, 
                               draw.lines = TRUE, 
                               lines.width = NULL, 
                               group.bar.height = 0.02, 
                               combine = TRUE) 
{
  cells <- cells %||% colnames(x = object)
  if (is.numeric(x = cells)) {
    cells <- colnames(x = object)[cells]
  }
  assay <- assay %||% DefaultAssay(object = object)
  DefaultAssay(object = object) <- assay
  features <- features %||% VariableFeatures(object = object)
  ## Why reverse???
  features <- rev(x = unique(x = features))
  disp.max <- disp.max %||% ifelse(test = slot == "scale.data", 
                                   yes = 2.5, no = 6)
  possible.features <- rownames(x = GetAssayData(object = object, 
                                                 slot = slot))
  if (any(!features %in% possible.features)) {
    bad.features <- features[!features %in% possible.features]
    features <- features[features %in% possible.features]
    if (length(x = features) == 0) {
      stop("No requested features found in the ", slot, 
           " slot for the ", assay, " assay.")
    }
    warning("The following features were omitted as they were not found in the ", 
            slot, " slot for the ", assay, " assay: ", paste(bad.features, 
                                                             collapse = ", "))
  }
  data <- as.data.frame(x = as.matrix(x = t(x = GetAssayData(object = object, 
                                                             slot = slot)[features, cells, drop = FALSE])))

  object <- suppressMessages(expr = StashIdent(object = object, 
                                               save.name = "ident"))
  group.by <- group.by %||% "ident"
  groups.use <- object[[c(group.by, additional.group.by)]][cells, , drop = FALSE]
  plots <- list()
  for (i in group.by) {
    data.group <- data
    group.use <- groups.use[, c(i, additional.group.by), drop = FALSE]
    
    for(colname in colnames(group.use)){
      if (!is.factor(x = group.use[[colname]])) {
        group.use[[colname]] <- factor(x = group.use[[colname]])
      }  
    }
    
    if (draw.lines) {
      lines.width <- lines.width %||% ceiling(x = nrow(x = data.group) * 
                                                0.0025)
      placeholder.cells <- sapply(X = 1:(length(x = levels(x = group.use[[i]])) * 
                                           lines.width), FUN = function(x) {
                                             return(Seurat:::RandomName(length = 20))
                                           })
      placeholder.groups <- data.frame(foo=rep(x = levels(x = group.use[[i]]), times = lines.width))
      placeholder.groups[additional.group.by] = NA
      colnames(placeholder.groups) <- colnames(group.use)
      rownames(placeholder.groups) <- placeholder.cells
      
      group.levels <- levels(x = group.use[[i]])

      group.use <- sapply(group.use, as.vector)
      rownames(x = group.use) <- cells
      
      group.use <- rbind(group.use, placeholder.groups)
      
      na.data.group <- matrix(data = NA, nrow = length(x = placeholder.cells), 
                              ncol = ncol(x = data.group), dimnames = list(placeholder.cells, 
                                                                           colnames(x = data.group)))
      data.group <- rbind(data.group, na.data.group)
    }
    
    #group.use = group.use[order(group.use[[i]]), , drop=F]
    group.use <- group.use[with(group.use, eval(parse(text=paste('order(', paste(c(i, additional.group.by), collapse=', '), ')', sep='')))), , drop=F]
    
    plot <- Seurat:::SingleRasterMap(data = data.group, raster = raster, 
                            disp.min = disp.min, disp.max = disp.max, feature.order = features, 
                            cell.order = rownames(x = group.use), group.by = group.use[[i]])
    
    if (group.bar) {
      pbuild <- ggplot_build(plot = plot)
      group.use2 <- group.use
      cols <- list()
      na.group <- Seurat:::RandomName(length = 20)
      for (colname in rev(x = colnames(group.use2))){
        if (colname == group.by){
          colid = paste0('Identity (', colname, ')')
        } else {
          colid = colname
        }
        
        if (draw.lines) {
          levels(x = group.use2[[colname]]) <- c(levels(x = group.use2[[colname]]), na.group)  
          group.use2[placeholder.cells, colname] <- na.group
          cols[[colname]] <- c(scales::hue_pal()(length(x = levels(x = group.use[[colname]]))), "#FFFFFF")
        } else {
          cols[[colname]] <- c(scales::hue_pal()(length(x = levels(x = group.use[[colname]]))))
        }
        names(x = cols[[colname]]) <- levels(x = group.use2[[colname]])
        
        
        y.range <- diff(x = pbuild$layout$panel_params[[1]]$y.range)
        y.pos <- max(pbuild$layout$panel_params[[1]]$y.range) + y.range * 0.015
        y.max <- y.pos + group.bar.height * y.range
        pbuild$layout$panel_params[[1]]$y.range <- c(pbuild$layout$panel_params[[1]]$y.range[1], y.max)
        
        plot <- suppressMessages(plot + 
                                annotation_raster(raster = t(x = cols[[colname]][group.use2[[colname]]]),  xmin = -Inf, xmax = Inf, ymin = y.pos, ymax = y.max) + 
                                annotation_custom(grob = grid::textGrob(label = colid, hjust = 0, gp = grid::gpar(cex = 0.75)), ymin = mean(c(y.pos, y.max)), ymax = mean(c(y.pos, y.max)), xmin = Inf, xmax = Inf) +
                                coord_cartesian(ylim = c(0, y.max), clip = "off")) 
        
        #temp <- as.data.frame(cols[[colname]][levels(group.use[[colname]])])
        #colnames(temp) <- 'color'
        #temp$x <- temp$y <- 1
        #temp[['name']] <- as.factor(rownames(temp))
        
        #temp <- ggplot(temp, aes(x=x, y=y, fill=name)) + geom_point(shape=21, size=5) + labs(fill=colname) + theme(legend.position = "bottom")
        #legend <- get_legend(temp)
        #multiplot(plot, legend, heights=3,1)
        
        if ((colname == group.by) && label) {
          x.max <- max(pbuild$layout$panel_params[[1]]$x.range)
          x.divs <- pbuild$layout$panel_params[[1]]$x.major
          group.use$x <- x.divs
          label.x.pos <- tapply(X = group.use$x, INDEX = group.use[[colname]],
                                FUN = median) * x.max
          label.x.pos <- data.frame(group = names(x = label.x.pos), 
                                    label.x.pos)
          plot <- plot + geom_text(stat = "identity", 
                                   data = label.x.pos, aes_string(label = "group", 
                                                                  x = "label.x.pos"), y = y.max + y.max * 
                                     0.03 * 0.5, angle = angle, hjust = hjust, 
                                   size = size)
          plot <- suppressMessages(plot + coord_cartesian(ylim = c(0, 
                                                                   y.max + y.max * 0.002 * max(nchar(x = levels(x = group.use[[colname]]))) * 
                                                                     size), clip = "off"))
        }
      }
    }
    plot <- plot + theme(line = element_blank())
    plots[[i]] <- plot
  }
  if (combine) {
    plots <- CombinePlots(plots = plots)
  }
  return(plots)
}

Hello,

This is not working for me. For the code with the looping, did we need to change anything? I left everything untouched in the long code and then just changed the one line that you provided first to reflect my object.

I keep getting this error:
Error in [.data.frame(groups.use, , c(i, additional.group.by), drop = FALSE) :
undefined columns selected

Thanks in advance to whomever helps me!

yanranz · 2021-03-08T22:57:15Z

Error: Must request at least one colour from a hue palette.
Called from: (scales::hue_pal())(length(x = levels(x = group.use[[colname]])))
Browse[1]> 1

Error in ((scales::hue_pal())(4))(length(x = levels(x = group.use[[colname]]))) :
attempt to apply non-function

elliefewings · 2021-03-16T10:28:58Z

Hi everyone,

I fixed the issues in the script above.

Looks like a few things were deprecated like pbuild$layout$panel_params[[1]]$x.major and there was some issues where levels were being called on the wrong object. It works now for me. I have posted it on my github. Obviously full credit goes to @arkal.

Hope that helps!

ondp15 · 2021-04-07T20:49:36Z

Hello,
Thanks for compiling arkal's code into a package. With the latest updates, I no longer can use the code and I'm trying your package. When I run the following command,

mesenchymalmarkers <- FindAllMarkers(mesenchymal, only.pos = TRUE, min.pct = 0.1, logfc.threshold = 0.10)
mesenchymalmarkers %>% group_by(cluster) %>% top_n(n = 75, wt = avg_log2FC)

Heatmapsample<- DoMultiBarHeatmap(mesenchymal, features = mesenchymalmarkers$gene, group.by = "CellType", additional.group.by = "samples" + theme(axis.text.y = element_text(size = 1.2)))
it results in the following warning:

Warning: CombinePlots is being deprecated. Plots should now be combined using the patchwork system.

I also tried removing the argument features and didn't work either. I also tried using the code above source code, which worked before and now it doesn't work either.

In the past year around March 2020, I used the function as described above and now it stopped working. I downloaded and updated the latest R version for mac as well as the packages.

Thanks for your help!

IceBear-nie · 2021-04-15T06:59:27Z

@ondp15 I had a similar issue -- my error was "could not find function CombinePlots". I just deleted the part about the CombinePlots, i.e.,
if (combine) {
plots <- CombinePlots(plots = plots)
}
and it allowed me to plot the right heatmap I wanted. I'm a newbie with R and seurat so I do not understand logically how deleting the portion affects the function, but at least the each gene/cell seems to show the right expression level still.

ktessema · 2021-04-15T07:06:57Z

Hello,
Thanks for compiling arkal's code into a package. With the latest updates, I no longer can use the code and I'm trying your package. When I run the following command,

mesenchymalmarkers <- FindAllMarkers(mesenchymal, only.pos = TRUE, min.pct = 0.1, logfc.threshold = 0.10)
mesenchymalmarkers %>% group_by(cluster) %>% top_n(n = 75, wt = avg_log2FC)

Heatmapsample<- DoMultiBarHeatmap(mesenchymal, features = mesenchymalmarkers$gene, group.by = "CellType", additional.group.by = "samples" + theme(axis.text.y = element_text(size = 1.2)))
it results in the following warning:

Warning: CombinePlots is being deprecated. Plots should now be combined using the patchwork system.

I also tried removing the argument features and didn't work either. I also tried using the code above source code, which worked before and now it doesn't work either.

In the past year around March 2020, I used the function as described above and now it stopped working. I downloaded and updated the latest R version for mac as well as the packages.

Thanks for your help!

Hey! I think you need to close the heatmap function parentheses before changing the theme. So add a ) after "samples" and remove one from the end of your code. There may be additional issues after fixing that though :P

elliefewings · 2021-04-15T07:31:03Z

Hi @iheartfoosball , @ondp15 also commented this on my own github page and I responded there. The CombinePlots deprecation message is only a warning. It does not (in my case) stop the package running and right now I do not have the time to change the plot combination function. If you have any other error message that is preventing the package from working please let me know!

ktessema · 2021-04-15T07:46:40Z

Hi @elliefewings - thanks for clarifying! I have not attempted this recently, just wanted to suggest that syntax edit for @ondp15. If the function has no other issues, moving the parenthesis should do the trick :)

xmc811 · 2021-04-18T17:01:29Z

I would like to recommend the R package Scillus I wrote for this functionality. It takes advantage of the ComplexHeatmap. The full vignette is at https://scillus.netlify.app/

Following are the examples:

plot_heatmap(dataset = scRNA_int, 
              markers = markers,
              sort_var = c("seurat_clusters","sample"),
              anno_var = c("seurat_clusters","sample","percent.mt","S.Score","G2M.Score"),
              anno_colors = list("Set2",                                             # RColorBrewer palette
                                 c("red","orange","yellow","purple","blue","green"), # color vector
                                 "Reds",
                                 c("blue","white","red"),                            # Three-color gradient
                                 "Greens"))

plot_heatmap(dataset = scRNA_int,
             n = 6,
             markers = markers,
             sort_var = c("seurat_clusters","sample"),
             anno_var = c("seurat_clusters","sample","percent.mt"),
             anno_colors = list("Set2",
                                c("red","orange","yellow","purple","blue","green"),
                                "Reds"),
             hm_limit = c(-1,0,1),
             hm_colors = c("purple","black","yellow"))

Sophia409 · 2021-12-23T12:30:48Z

I also get the following error and have no idea how to fix it. How did you solve it?

DoMultiBarHeatmap(PVN.neuron,assay = "RNA",features = gene, group.by='id',additional.group.by = 'batch')
 Error: Must request at least one colour from a hue palette.

@iheartfoosball @arkal @xmc811 @elliefewings @MinBio

martin-jeremy · 2022-02-03T13:41:45Z

I also get the following error and have no idea how to fix it. How did you solve it?
DoMultiBarHeatmap(PVN.neuron,assay = "RNA",features = gene, group.by='id',additional.group.by = 'batch')
 Error: Must request at least one colour from a hue palette.
@iheartfoosball @arkal @xmc811 @elliefewings @MinBio

same, I can't use this function and I don't understand why.
It block at line 110 cols[[colname]] <- c(scales::hue_pal()(length(x = levels(x = group.use[[colname]]))), "#FFFFFF") because length(x = levels(x = group.use[[colname)) = 0 so that is calling sclaes::hue_pal()(0) ...

I try to pass over it with draw.lines = F, but it's the same problem line 112...

Any suggestion ?

JCMolBio · 2022-02-11T16:16:23Z

I also get the following error and have no idea how to fix it. How did you solve it?
DoMultiBarHeatmap(PVN.neuron,assay = "RNA",features = gene, group.by='id',additional.group.by = 'batch')
 Error: Must request at least one colour from a hue palette.
@iheartfoosball @arkal @xmc811 @elliefewings @MinBio
same, I can't use this function and I don't understand why. It block at line 110 cols[[colname]] <- c(scales::hue_pal()(length(x = levels(x = group.use[[colname]]))), "#FFFFFF") because length(x = levels(x = group.use[[colname)) = 0 so that is calling sclaes::hue_pal()(0) ...

I try to pass over it with draw.lines = F, but it's the same problem line 112...

Any suggestion ?

I found that adding the argument "group.bar=T" at least gets me to display the heatmap.

Pengwei-Ren · 2022-07-09T20:33:12Z

DoMultiBarHeatmap is easy to use for me who don't have coding experience. It did give more details on group bar.

I tried to set specific color for the group bars for the seurat clusters and conditions using script below:

cols.use <- list(
seurat_clusters = c('0' = "#F8766D", '1' = "#7CAE00", '2' = "#00BFC4", '3' = "#C77CFF"),
orig.ident = c(KO1 = "#F8766D", KO2 = "#7CAE00", WT1 = "#00BFC4", WT2 = "#C77CFF")
)
DoMultiBarHeatmap(seuratobject, assay = 'SCT',
features = features,
group.by= 'seurat_clusters',
additional.group.by = 'orig.ident',
cols.use = cols.use,
size = 3)

While R reports "unused argument (cols.use = cols.use)", I also tried "group.colors" also not working.

Then the color is the same as without adding "cols.use" (cluster color is the default one of UMAP, while sample color is new, wanted to change it as in the UMAP).

Anyone can help? Thank you!

JonathanNoonan · 2022-07-22T00:16:55Z

Hi All - This works perfectly with single cell heatmaps, but the grouping bars end up offset if you try to plot averaged expression. Does anyone have ant idea how to fix this? - Definitely wish this was added to Seurat - great widget.

Elo-mars · 2022-11-10T08:47:33Z

Thanks for the code @arkal . It's very helpful. I noticed that the legend is missing for the new group bars. Could you add it to the heatmap? thank you.

Indeed, this is what I was looking for too, did you manage it?

thanks for creating this @arkal ! this is very useful for people in my group
but at the moment, I have to go back to the metadata to see how many cells of each secondary group I have, so I know which color or the second bar is corresponding to what (the additional group)

did you figure out how to add a legend for this additional group?

Thanks

Elodie

TauHHCC · 2023-03-31T17:36:10Z

Hi @arkal ,
Thank you, that would be really helpful.
I'm trying to use your updated function to run, and I got this error

DoMultiBarHeatmap(hcf.ab.all,features = top25$gene, group.by='seurat_clusters', additional.group.by = 'Timepoints')
Error in gpar(cex = 0.75) : could not find function "gpar"
Called from: grob(label = label, x = x, y = y, just = just, hjust = hjust, 
    vjust = vjust, rot = rot, check.overlap = check.overlap, 
    name = name, gp = gp, vp = vp, cl = "text")

I also tried change additional.group.by = to additional.group.sort.by = and got this errors

DoMultiBarHeatmap(hcf.ab.all,features = top25$gene, group.by='seurat_clusters', additional.group.sort.by = 'Timepoints')
Error during wrapup: could not find function "gpar"
Error: no more error handlers available (recursive errors?); invoking 'abort' restart

Could you please help to point_

Best

liviu- · 2023-09-27T13:52:42Z

draw.lines = F and label = F fixed for me the errors for the DoMultiBarHeatmap() function above.

denvercal1234GitHub · 2023-11-15T11:11:06Z

I would like to recommend the R package Scillus I wrote for this functionality. It takes advantage of the ComplexHeatmap. The full vignette is at https://scillus.netlify.app/

Following are the examples:

plot_heatmap(dataset = scRNA_int, 
              markers = markers,
              sort_var = c("seurat_clusters","sample"),
              anno_var = c("seurat_clusters","sample","percent.mt","S.Score","G2M.Score"),
              anno_colors = list("Set2",                                             # RColorBrewer palette
                                 c("red","orange","yellow","purple","blue","green"), # color vector
                                 "Reds",
                                 c("blue","white","red"),                            # Three-color gradient
                                 "Greens"))

plot_heatmap(dataset = scRNA_int,
             n = 6,
             markers = markers,
             sort_var = c("seurat_clusters","sample"),
             anno_var = c("seurat_clusters","sample","percent.mt"),
             anno_colors = list("Set2",
                                c("red","orange","yellow","purple","blue","green"),
                                "Reds"),
             hm_limit = c(-1,0,1),
             hm_colors = c("purple","black","yellow"))

Hi @xmc811 --- Does this heatmap group the cells by hierarchical clustering for a metadata variable? Or is there a way to do that here?

Thank you so much

kvn95ss · 2023-12-11T09:42:28Z

Hello!
With Seurat v5, the scripts don't seem to work any more. I tried to fix it but I wasn't successful :(

If anyone can have a go at fixing the script it would be great! For now I'll use Seurat v4 to generate the heatmap.

nghiaplt1611 · 2024-04-16T16:42:51Z

FWIW, I modified DoHeatmap to do this in a simple way. I didn't get to adding a legend for the additional groups but the bar on the top can be seen to better understand your data.

E.g.: MySeuratObj integrated 3 samples together and i wanted to see the DGE per identified cluster (seurat_idents) and additionally wanted to see the sample each column in the heatmap came from (orig.ident) so I can do

> DoMultiBarHeatmap(MySeuratObj, assay = 'SCT', features = top_10$gene, group.by='seurat_idents', additional.group.by = 'orig.ident')

This is the result (genes and clusters, etc are redacted since this is unpublished work ;) )

Here's the code

suppressPackageStartupMessages({
  library(rlang)
})

DoMultiBarHeatmap <- function (object, 
                               features = NULL, 
                               cells = NULL, 
                               group.by = "ident", 
                               additional.group.by = NULL, 
                               group.bar = TRUE, 
                               disp.min = -2.5, 
                               disp.max = NULL, 
                               slot = "scale.data", 
                               assay = NULL, 
                               label = TRUE, 
                               size = 5.5, 
                               hjust = 0, 
                               angle = 45, 
                               raster = TRUE, 
                               draw.lines = TRUE, 
                               lines.width = NULL, 
                               group.bar.height = 0.02, 
                               combine = TRUE) 
{
  cells <- cells %||% colnames(x = object)
  if (is.numeric(x = cells)) {
    cells <- colnames(x = object)[cells]
  }
  assay <- assay %||% DefaultAssay(object = object)
  DefaultAssay(object = object) <- assay
  features <- features %||% VariableFeatures(object = object)
  ## Why reverse???
  features <- rev(x = unique(x = features))
  disp.max <- disp.max %||% ifelse(test = slot == "scale.data", 
                                   yes = 2.5, no = 6)
  possible.features <- rownames(x = GetAssayData(object = object, 
                                                 slot = slot))
  if (any(!features %in% possible.features)) {
    bad.features <- features[!features %in% possible.features]
    features <- features[features %in% possible.features]
    if (length(x = features) == 0) {
      stop("No requested features found in the ", slot, 
           " slot for the ", assay, " assay.")
    }
    warning("The following features were omitted as they were not found in the ", 
            slot, " slot for the ", assay, " assay: ", paste(bad.features, 
                                                             collapse = ", "))
  }
  data <- as.data.frame(x = as.matrix(x = t(x = GetAssayData(object = object, 
                                                             slot = slot)[features, cells, drop = FALSE])))

  object <- suppressMessages(expr = StashIdent(object = object, 
                                               save.name = "ident"))
  group.by <- group.by %||% "ident"
  groups.use <- object[[c(group.by, additional.group.by)]][cells, , drop = FALSE]
  plots <- list()
  for (i in group.by) {
    data.group <- data
    group.use <- groups.use[, c(i, additional.group.by), drop = FALSE]
    
    for(colname in colnames(group.use)){
      if (!is.factor(x = group.use[[colname]])) {
        group.use[[colname]] <- factor(x = group.use[[colname]])
      }  
    }
    
    if (draw.lines) {
      lines.width <- lines.width %||% ceiling(x = nrow(x = data.group) * 
                                                0.0025)
      placeholder.cells <- sapply(X = 1:(length(x = levels(x = group.use[[i]])) * 
                                           lines.width), FUN = function(x) {
                                             return(Seurat:::RandomName(length = 20))
                                           })
      placeholder.groups <- data.frame(foo=rep(x = levels(x = group.use[[i]]), times = lines.width))
      placeholder.groups[additional.group.by] = NA
      colnames(placeholder.groups) <- colnames(group.use)
      rownames(placeholder.groups) <- placeholder.cells
      
      group.levels <- levels(x = group.use[[i]])

      group.use <- sapply(group.use, as.vector)
      rownames(x = group.use) <- cells
      
      group.use <- rbind(group.use, placeholder.groups)
      
      na.data.group <- matrix(data = NA, nrow = length(x = placeholder.cells), 
                              ncol = ncol(x = data.group), dimnames = list(placeholder.cells, 
                                                                           colnames(x = data.group)))
      data.group <- rbind(data.group, na.data.group)
    }
    
    #group.use = group.use[order(group.use[[i]]), , drop=F]
    group.use <- group.use[with(group.use, eval(parse(text=paste('order(', paste(c(i, additional.group.by), collapse=', '), ')', sep='')))), , drop=F]
    
    plot <- Seurat:::SingleRasterMap(data = data.group, raster = raster, 
                            disp.min = disp.min, disp.max = disp.max, feature.order = features, 
                            cell.order = rownames(x = group.use), group.by = group.use[[i]])
    
    if (group.bar) {
      pbuild <- ggplot_build(plot = plot)
      group.use2 <- group.use
      cols <- list()
      na.group <- Seurat:::RandomName(length = 20)
      for (colname in rev(x = colnames(group.use2))){
        if (colname == group.by){
          colid = paste0('Identity (', colname, ')')
        } else {
          colid = colname
        }
        
        if (draw.lines) {
          levels(x = group.use2[[colname]]) <- c(levels(x = group.use2[[colname]]), na.group)  
          group.use2[placeholder.cells, colname] <- na.group
          cols[[colname]] <- c(scales::hue_pal()(length(x = levels(x = group.use[[colname]]))), "#FFFFFF")
        } else {
          cols[[colname]] <- c(scales::hue_pal()(length(x = levels(x = group.use[[colname]]))))
        }
        names(x = cols[[colname]]) <- levels(x = group.use2[[colname]])
        
        
        y.range <- diff(x = pbuild$layout$panel_params[[1]]$y.range)
        y.pos <- max(pbuild$layout$panel_params[[1]]$y.range) + y.range * 0.015
        y.max <- y.pos + group.bar.height * y.range
        pbuild$layout$panel_params[[1]]$y.range <- c(pbuild$layout$panel_params[[1]]$y.range[1], y.max)
        
        plot <- suppressMessages(plot + 
                                annotation_raster(raster = t(x = cols[[colname]][group.use2[[colname]]]),  xmin = -Inf, xmax = Inf, ymin = y.pos, ymax = y.max) + 
                                annotation_custom(grob = grid::textGrob(label = colid, hjust = 0, gp = grid::gpar(cex = 0.75)), ymin = mean(c(y.pos, y.max)), ymax = mean(c(y.pos, y.max)), xmin = Inf, xmax = Inf) +
                                coord_cartesian(ylim = c(0, y.max), clip = "off")) 
        
        #temp <- as.data.frame(cols[[colname]][levels(group.use[[colname]])])
        #colnames(temp) <- 'color'
        #temp$x <- temp$y <- 1
        #temp[['name']] <- as.factor(rownames(temp))
        
        #temp <- ggplot(temp, aes(x=x, y=y, fill=name)) + geom_point(shape=21, size=5) + labs(fill=colname) + theme(legend.position = "bottom")
        #legend <- get_legend(temp)
        #multiplot(plot, legend, heights=3,1)
        
        if ((colname == group.by) && label) {
          x.max <- max(pbuild$layout$panel_params[[1]]$x.range)
          x.divs <- pbuild$layout$panel_params[[1]]$x.major
          group.use$x <- x.divs
          label.x.pos <- tapply(X = group.use$x, INDEX = group.use[[colname]],
                                FUN = median) * x.max
          label.x.pos <- data.frame(group = names(x = label.x.pos), 
                                    label.x.pos)
          plot <- plot + geom_text(stat = "identity", 
                                   data = label.x.pos, aes_string(label = "group", 
                                                                  x = "label.x.pos"), y = y.max + y.max * 
                                     0.03 * 0.5, angle = angle, hjust = hjust, 
                                   size = size)
          plot <- suppressMessages(plot + coord_cartesian(ylim = c(0, 
                                                                   y.max + y.max * 0.002 * max(nchar(x = levels(x = group.use[[colname]]))) * 
                                                                     size), clip = "off"))
        }
      }
    }
    plot <- plot + theme(line = element_blank())
    plots[[i]] <- plot
  }
  if (combine) {
    plots <- CombinePlots(plots = plots)
  }
  return(plots)
}

Hello, I have run this code and it worked. But I need to know more about the data that are grouping in the additional.group.by, I want to show these data the same as the way that group.by (Identity) is showing and plotting (additional labelling for columns used in additional.group.by). How can I perform this? Thanks a lot.

satijalab closed this as completed Oct 27, 2019

Dragonmasterx87 mentioned this issue Sep 28, 2020

label-annotated heatmap using DoHeatmap #3554

Closed

gabrieljm00 mentioned this issue Apr 2, 2021

DoHeatmap clustering specific genes and not top x genes #2261

Closed

ChristinaSteyn mentioned this issue May 16, 2022

Which assay and slot to use for visualising gene expression in different scenarios #5966

Closed

Multiple group bars on heatmap #2201

Multiple group bars on heatmap #2201

Comments

zehualilab commented Oct 14, 2019

satijalab commented Oct 27, 2019

arkal commented Nov 14, 2019 • edited Loading

HomairaH commented Dec 13, 2019 • edited Loading

arkal commented Dec 17, 2019

HomairaH commented Dec 29, 2019

MinBio commented Jan 15, 2020 • edited Loading

arkal commented Jan 16, 2020

kaizen89 commented Mar 26, 2020

ktessema commented Mar 27, 2020

kaizen89 commented Mar 27, 2020

csung331 commented Jul 16, 2020

ktessema commented Jul 16, 2020

csung331 commented Jul 16, 2020

ktessema commented Jul 16, 2020

csung331 commented Jul 16, 2020

ktessema commented Jul 16, 2020 • edited Loading

csung331 commented Jul 16, 2020

ktessema commented Jul 16, 2020

PrashINRA commented Jul 29, 2020

Diozao666 commented Aug 2, 2020 • edited Loading

PrashINRA commented Aug 2, 2020

shaln commented Feb 2, 2021

ktessema commented Feb 4, 2021

InMyGenes commented Feb 11, 2021

yanranz commented Mar 8, 2021 • edited Loading

elliefewings commented Mar 16, 2021

ondp15 commented Apr 7, 2021

IceBear-nie commented Apr 15, 2021

ktessema commented Apr 15, 2021

elliefewings commented Apr 15, 2021

ktessema commented Apr 15, 2021

xmc811 commented Apr 18, 2021

Sophia409 commented Dec 23, 2021 • edited Loading

martin-jeremy commented Feb 3, 2022

JCMolBio commented Feb 11, 2022

Pengwei-Ren commented Jul 9, 2022

JonathanNoonan commented Jul 22, 2022

Elo-mars commented Nov 10, 2022

TauHHCC commented Mar 31, 2023

liviu- commented Sep 27, 2023

denvercal1234GitHub commented Nov 15, 2023

kvn95ss commented Dec 11, 2023

nghiaplt1611 commented Apr 16, 2024

arkal commented Nov 14, 2019 •

edited

Loading

HomairaH commented Dec 13, 2019 •

edited

Loading

MinBio commented Jan 15, 2020 •

edited

Loading

ktessema commented Jul 16, 2020 •

edited

Loading

Diozao666 commented Aug 2, 2020 •

edited

Loading

yanranz commented Mar 8, 2021 •

edited

Loading

Sophia409 commented Dec 23, 2021 •

edited

Loading