using heat map to measure gene expression at different developmental stages

Question

[Kai_Qi]

Hi All:

I have a matrix that contains the gene EntrezID and the counts number at different developmental stage, the head of which looks like this:

> head(count_cts)
  E11_rep1 E11_rep2 E14_rep1 E14_rep2 E18_rep1 E18_rep2 Adult_rep1 Adult_rep2
21749     4638     3254     1472     1862      632     1019        799        684
98403     8794     7109     5282     7475     3778     5853       2569       2595
18393     7917     7888     3758     5916     2814     4480       2997       2752
21961     6972     9219     6400     5051     2583     4079       5272       5210
76709    15345    19298    16327    21389    19766    28117      14388      13096
67534    12314     9867     5107     7963     4060     6005       1404       1397

It contains the developmental stages, each contains 2 replicates, the rowname of the the matrix is the geneID. The number is the counts I got from Rsubread. The geneID are already selected using a criteria. Now I want to use heatmap to describe the gene expression level change.

My question is: How should I process the counts before feed them into heatmap?

Thank you,

[Kai_Qi]

I tried. It gave a heatmap with all black color.

[N15]

You could try adding another break in. I would keep it simple and call the heatmap on your normalized count matrix. After that, you can modify the parameters one by one to change colors and scale bar breaks to your liking.

library(pheatmap)
pheatmap(counts)

[Kai_Qi]

Will try. Your advice looks working better. Thanks!

[Kai_Qi]

Hi: I have read some short manuals for heatmap for getting better images but I failed to get one with color showing the dynamic change of reads counts(normalized).

my_palette <- colorRampPalette(c("red", "yellow", "green"))(n = 299)

# (optional) defines the color breaks manually for a "skewed" color transition
col_breaks = c(seq(25,500,length=100),  # for red
               seq(20.01,25,length=150),           # for yellow
               seq(0,20,length=100))             # for green

# creates a 5 x 5 inch image
png("heatmaps_in_r.png",    # create PNG for the heat map        
    width = 5*300,        # 5 x 300 pixels
    height = 5*300,
    res = 300,            # 300 pixels per inch
    pointsize = 8)        # smaller font size
heatmap.2(count_cts, cellnote = count_cts, main = "Dynamic_Expression", density.info="none", trace="none",margins =c(12,20), col=my_palette,  dendrogram = "row", Colv="NA")

The heatmap generated is showing with single color. I guess it is because the reads counts value show too big variation for different genes. I am wondering how can I set the color to show the change at different stages:

> head(count_cts)
       E11_rep1  E11_rep2  E14_rep1  E14_rep2   E18_rep1  E18_rep2 Adult_rep1 Adult_rep2
21749  36.30160  25.52971  12.69430  12.19242   6.178822   6.58037   9.699078   8.476963
98403  68.83059  55.77466  45.55116  48.94647  36.936061  37.79677  31.185147  32.160409
18393  61.96632  61.88641  32.40842  38.73810  27.511402  28.93038  36.380648  34.106145
21961  54.56981  72.32896  55.19262  33.07406  25.253003  26.34085  63.996922  64.568683
76709 120.10523 151.40517 140.80155 140.05567 193.244623 181.57044 174.656244 162.301625
67534  96.38161  77.41293  44.04199  52.14191  39.693067  38.77834  17.043186  17.313330

For example for gene 21749, 36.30160 is high at E11 compared to 8.476963 at Adult; however, this value is the low expression in gene 98403. How to solve this conflicts? Also after heatmap, there will be clustering. Is there any methods to get the information on which genes are in the same cluster?

[N15]

Try log transforming your count matrix before plotting with heatmap. You can export log-transformed counts from edgeR directly:

tmm<-cpm(dge, normalized.lib.size=T, log=T)

[Kai_Qi]

Thanks for the reply. Though the numbers looks better, but it still is a single colored heatmap.

How can I show you the map I got to make my expression better to be understood.

Thanks a lot,

Answer 1

[N15]

I would normalize counts using either deseq2 or edgeR, and then plot them in using a heatmap tool that can cluster columns based on similarity. That should show you major differences in EntrezID/function between samples.

https://genviz.org/module-04-expression/0004/02/01/DifferentialExpression/

[Kai_Qi]

I know edgeR has a function:

x <- calNormFactors(x, method="TMM")

I will get the normalized factor for each group. I am wondering should I multiply each reads number with their corresponding normalize factors?

Thank you,

[N15]

No let edgeR do it for you!!!!

dge<-DGEList(counts,group = c(1,2,3,4)) 

dge<-calcNormFactors(dge)

tmm<-cpm(dge, normalized.lib.size=T, log=F) #non-LOG TMM normalized counts

[Kai_Qi]

I almost got there but got an error at the very last step:

my_palette <- colorRampPalette(c("red", "yellow", "blue"))(n = 299)

col_breaks = c(seq(0,20,length=100),  # for red
               seq(20.01,50,length=150),           # for yellow
               seq(50.01,300,length=100))             # for blue

png("heatmaps_in_r.png",    # create PNG for the heat map        
    width = 5*300,        # 5 x 300 pixels
    height = 5*300,
    res = 300,            # 300 pixels per inch
    pointsize = 8)

every thing looked fine, until the following step:

heatmap.2(mat_data, cellnote = count_cts, main = "Dynamic_Expression", notecol = "black", density.info="none", trace="none",margins =c(12,20), col=my_palette, breaks = col_breaks, dendrogram = "row", Colv="NA")

Error in image.default(1:nc, 1:nr, x, xlim = 0.5 + c(0, nc), ylim = 0.5 + : must have one more break than colour

what can be the mistakes, is that because the value of the normalized counts exceed the color breaks?

[N15]

Perhaps. Did you try it without col_breaks?