This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How to scale raw Counts to obtain z-score Heatmap in R?

Hi everybody, im trying to replicate this image from a paper (of which I have replicated all analysis except for this heatmap) :

enter image description here

Caption from the paper : Heatmap of unsupervised clustering of platelet mRNA profiles of healthy donors (red) and patients with cancer (gray). Paper link : https://www.cell.com/cancer-cell/fulltext/S1535-6108(15)00349-9

What I did is the following :

data <- cpm(counts,prior.count = 1,log=TRUE)
data.scaled <- scale(data,center = TRUE, scale = TRUE) #I tried all different scale types
data.scaled <- data.scaled[rownames(TopGenes),] #Choose just some genes so they are not all as the paper did
h <- pheatmap(normList.scaled[0:300,], cluster_rows=FALSE,name ="Z-score",
                cluster_cols=FALSE,annotation_col=info)#,annotation_col=info) #0:300 because all of them make my computer slow

And this is what comesenter image description here out :

Now the order of the samples is not important ( I did another order type ) , but i don't understand how they computed the z-scores. They didnt specify the steps so i guessed they first did the cpm and then scaled? But it doesn't look right.

geo r
• 5,329 views • link
written by ali • ali

That looks like a ComplexHeatmap with a topAnnotation. CH clusters using hclust(dist(mat)) by default. You may want to try using CH instead of pheatmap.

• 0 views • link
written by Ram Ram

I'd recommend emailing the authors and requesting the exact code used to generate the heatmap.

• 0 views • link
written by Ram Ram

Thank you for the reply, I tried row scaling too but still is completely different from their work. I have seen the supplemental files but is not said how to do it still. They just show another Heatmap saying that refers to the count per millions but I did it too (cpm and then scaling) but my values are way higher than 1.5 (like 10)

• 0 views • link
written by ali • ali

1 answer

See my answer here on how-to and why one scales data for heatmaps Scaling RNA-Seq data before clustering?

• 0 views • link
written by ATpoint ATpoint

Hey thanks for this link, it was really useful. The preprocessing step is the same I did but with filtering FDR < 0.001 because thats how they did it and logCPM > 3. I have retrived something like this :

enter image description here

It looks really similar now (just inverted , maybe the gene order is exactly my opposite idk). Could you explain to me why I do have 4 and -4 and they have 1.5 and -1.5? I read what you wrote on the other post and I got what z-score is but , if the data I'm using is the same , how is possible that the z-scores differ. Thats the code I use now (still have to apply fisher even if i dont know how) :

fit <- glmQLFit(list,design,robust = TRUE) #list is the DGEList
single.contrast <- makeContrasts(Tumor  - HC, levels = design)
test.anov <- glmQLFTest(fit,contrast = single.contrast)



Filtering data by logCPM > 3 and FDR < 0.001 #####

tagsTblANOVA <- topTags(test.anov, n=nrow(test.anov$table))$table
tagsTblANOVA.logcpm <- tagsTblANOVA[tagsTblANOVA$logCPM>3,]
tagsTblAnova.fdr <- tagsTblANOVA.logcpm[tagsTblANOVA.logcpm$FDR <0.001,]

column_ha = HeatmapAnnotation(condition = info$Condition)
DEGs <- rownames(tagsTblAnova.fdr) #filtered by FDR after fit and test.
cts <- cpm(list, log = TRUE)[rownames(list) %in% DEGs,] #list is the DGELIST obj
jpeg(file="h1.jpg")
Heatmap(matrix = t(scale(t(cts))),
           cluster_rows = TRUE,
           clustering_method_rows = "ward.D2",
           cluster_columns = TRUE,
           clustering_method_columns = "ward.D2",
           show_row_names = FALSE, 
           show_column_names = FALSE,
           column_title = "logcounts scaled",
    top_annotation = column_ha)
dev.off()
• 0 views • link
written by ali • ali

Unless you do exactly the same code with exactly the same software versions and exactly the same everything, results will have nuance differences. Maybe they trimmed away some outliers from the heatmap color scale, can be many reasons. The message should be the same for the sake of reproducibility, that is what matters.

• 0 views • link
written by ATpoint ATpoint

thx before i close the post as solved, do you have any advise on how to implement fisher test on clustering into the Heatmap?

• 0 views • link
written by ali • ali

Log in to answer this question.

More posts like this

No similar posts found.