I am analysing RNA-Seq data for a set of genes from different human brain regions using data from the Allen Brain Atlas. I have made a heatmap using GraphPad Prism. In this heatmap, I show the average RPKM value for each gene of interest for each brain region (the average is taken as there are multiple samples for each brain region).

I am new to creating heatmaps but I know that you can use different methods to cluster the heatmap so that the data is more easy to interpret. I want to cluster the heatmap so that the brain regions which have the highest level of gene expression overall are grouped together.

I have seen online that R can also be used to create a heatmap, and I have had previous experience with R. Does anyone have any recommendations on any good tutorials (that are not too difficult) on creating heat-maps in R as well as clustering the data?


1 Answer 1


For visualization of expression values it usually makes sense to standardize the data first. The problem with plain expression values is that they differ a lot between genes and as such highly-expressed genes dominate the color scale. As you see in your heatmap, the first and last rows barely show any difference because the color scale is dominated by a single gene CHL1. Therefore the clustering would be due to expression level and not due to differences between columns. For this we standardize the data. Assuming your expression matrix had RPKMs (or similar values) we first log2-transform it and then standardize:

mat_log2 <- log2(mat+1)

#/ standardize (Z = Z-scored)
mat_Z <- t(scale(t(mat_log2)))

#/ cluster rows:
hclust_row <- hclust(dist(mat_Z), method="ward.D2")

#/ cluster columns (optional)
hclust_col <- hclust(dist(t(mat_Z)), method="ward.D2")

#/ plot heatmap:

hm <- 
          name="relative expression")

See the ComplexHeatmap documentation for details on the package itself, it is very powerful but also monolithic. This is just an example, any other heatmap package will do fine as well.

Important is to standardize data to put emphasis on differences between columns rather than on absolute expression levels.

  • $\begingroup$ Thanks for your reply. Would it be possible to know why the RPKM values are log2-transformed before the Z-score standardisation is done? $\endgroup$
    – ceno980
    Dec 28, 2021 at 9:10
  • $\begingroup$ @ceno980 log2 is generally done on genomics data to dampen the effects of outliers as the gene expression values have such a wide dynamic range. As the standardization divides internally by the standard deviation it just helps that single values skew the calculation overly much e.g. in case where some samples have 0 and the others have thousands of counts. The log2 scale tries to make things a bit more homoskedastic. I guess that is more important if you do things like PCA on the expression values directly rather than Z-scoring but for consistency I always do it. $\endgroup$
    – user3051
    Dec 28, 2021 at 9:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.