Add Heatmap

[ktrns]

It would be nice to add a tab for a heatmap. For now, this can be the Seurat Heatmap, but we might as well change the code to plot a nicer heatmap. (I have some code)

This is actually a common request at our facility, and I haven't thought about this before. We can also add this later, if it is outside the scope of Marius project.

Best & thanks Katrin

[mariusrueve]

I can create the tab and plots for the heatmap. Could you paste the code in this issue?

[ktrns]

This is the code we use in the main scrnaseq script:

# Heatmap of all differentially expressed genes
p = Seurat::DoHeatmap(sc, features=degs_filt$all$gene, group.colors=param$col_clusters) + 
  NoLegend() + 
  ggtitle("Heatmap of scaled gene expression data, all genes differentially expressed between a cluster and the rest")
p

My current client would like to have more freedom with regard to adding a dendrogram, changing colours and so on, so I wrote some code for him to do this in R:

# Extract raw, normalised and scaled data for all genes of interest
counts = GetAssayData(sc, assay="RNA", slot="counts")[genes,]
data = GetAssayData(sc, assay=param$normalisation_default, slot="data")[genes,]
scale.data = GetAssayData(sc, assay=param$normalisation_default, slot="scale.data")[genes,]

# Heatmap 
# Define colors for cells
cell_ident = sc[[]]$seurat_clusters
idents = sort(unique(cell_ident))
idents_col = ggsci::pal_npg()(length(idents))
names(idents_col) = idents
cell_col = idents_col[cell_ident]

# Define colors for gene expression values
col_fun = circlize::colorRamp2(c(-1.5,0,1.5), c("blue", "white", "yellow")) 

# Define annotation for cells 
top_annotation = ComplexHeatmap::HeatmapAnnotation(Identity=cell_ident, col=list(Identity=idents_col))

# Plot heatmap
ComplexHeatmap::Heatmap(matrix=scale.data, 
                        col=col_fun,
                        row_title="Genes",
                        column_title="Cells",
                        cluster_rows=FALSE,
                        cluster_columns=FALSE,
                        show_row_names=FALSE,
                        show_column_names=FALSE,
                        column_split=cell_ident, 
                        top_annotation=top_annotation,
                        heatmap_legend_param = list(title="Gene expression")
                        )

I suggest you start with the simpler Seurat code (first code chunk).

[mariusrueve]

@ktrns

plots_Heatmap <<- #global assignment
        Seurat::DoHeatmap(
          sc(), 
          features=unlist(strsplit(input$genes, "_"))[c(T, F)], #extract genes from selection
          group.colors=param$col_clusters) + 
        NoLegend()

produces following error: Warning: Error in Seurat::DoHeatmap: No requested features found in the scale.data slot for the SCT assay.

[mariusrueve]

@ktrns creating a heatmap works now. I used this code:

      plots_Heatmap <<-
        Seurat::DoHeatmap(
          sc(),
          features=unlist(strsplit(input$genes, "_"))[c(T, F)],
          group.colors=param$col_clusters,
          slot = "counts") +
        NoLegend()

which will result in this image:

Warning:

Warning in Seurat::DoHeatmap(sc(), features = unlist(strsplit(input$genes,  :
  The following features were omitted as they were not found in the counts slot for the SCT assay: LINC01342
Warning: Could not find LINC01342 in the default search locations, found in RNA assay instead

The gene is not shown in the plot. How can i fix this problem?

[ktrns]

Hi @mariusrueve,

Great that the heatmap is now part of the app. If I am not mistaken, the heatmap function uses the scale.data slot. The respective data set consists of variable genes only. It MAY be that this is different for the SCT and RNA assay, and we have changed this part a few times.

You can double-check this by dim(GetAssayData(sc, assay="SCT", slot="scale.data")), and dim(GetAssayData(sc, assay="RNA", slot="scale.data")). Can you try that and let me know the result? Also try the other slots ("counts" and "data").

Best and thanks Katrin

[mariusrueve]

These are the results @ktrns, I hope this is what you need.

> library(Seurat)
> sc = readRDS(file = "pbmc_2020-06-09.rds")

# scale.data
> dim(GetAssayData(sc, assay="SCT", slot="scale.data"))
[1] 3000 1113
> dim(GetAssayData(sc, assay="RNA", slot="scale.data"))
[1] 16220  1113

# counts
> dim(GetAssayData(sc, assay="SCT", slot="counts"))
[1] 13994  1113
> dim(GetAssayData(sc, assay="RNA", slot="counts"))
[1] 16220  1113

# data
> dim(GetAssayData(sc, assay="SCT", slot="data"))
[1] 13994  1113
> dim(GetAssayData(sc, assay="RNA", slot="data"))
[1] 16220  1113

[ktrns]

Hi @mariusrueve,

Yes, this is what we needed. Now look at the dimensions of scale.data for SCT and RNA. For SCT, we only kept 3,000 genes (variable genes). For RNA, we kept all genes.

This is something we solved in our current, soon to be submitted version of scrnaseq. So I'd think that you can safely ignore this warning.

However, you see there is still a difference between RNA and SCT, so in theory it might still happen that a gene is found in RNA but not in SCT. I don't think there is anything you can do about it except letting the user know.

Best wishes and a nice day Katrin