satijalab / seurat

R toolkit for single cell genomics
http://www.satijalab.org/seurat
Other
2.24k stars 902 forks source link

Multiple group bars on heatmap #2201

Closed zehualilab closed 4 years ago

zehualilab commented 4 years ago

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 commented 4 years ago

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 commented 4 years ago

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 ;) )

foobar

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 commented 4 years ago

@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:

1) 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?

2) 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 ;) )

foobar

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 commented 4 years ago

@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'))

image

# 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'))

image

# 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'))

image

# 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)

image

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)

image

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 commented 4 years ago

@arkal thanks so much again!

MinBio commented 4 years ago

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 commented 4 years ago

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 commented 4 years ago

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 commented 4 years ago

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 commented 4 years ago

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

csung331 commented 4 years ago

@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 commented 4 years ago

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

csung331 commented 4 years ago

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

ktessema commented 4 years ago

@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 commented 4 years ago

@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 commented 4 years ago

@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 commented 4 years ago

@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 commented 4 years ago

@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 commented 4 years ago

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 commented 4 years ago

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 commented 4 years ago

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 commented 3 years ago

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 commented 3 years ago

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 commented 3 years ago

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 ;) )

foobar

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 commented 3 years ago

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 commented 3 years ago

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 commented 3 years ago

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 commented 3 years ago

@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 commented 3 years ago

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 commented 3 years ago

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 commented 3 years ago

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 commented 3 years ago

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"))

image

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"))

image

Sophia409 commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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. Screenshot 2022-07-22 at 10 16 25

Elo-mars commented 1 year ago

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 commented 1 year ago

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- commented 11 months ago

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

denvercal1234GitHub commented 10 months ago

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"))

image

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"))

image

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 commented 9 months ago

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 commented 5 months ago

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 ;) )

foobar

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.