I want to make a heatmap by pheatmap package and I am new to this. I have a list of DEGs obtained by RNAseq analysis using DESeq2 and I only have normalized counts from DESeq2 and also FPKM values for these genes. To my knowledge, I need to plot Z-score values from the 3 replicates per sample for each gene (rows) but I dont know how to get Z-scores. Should I calculate the Z-scores first and then read it into R as a data frame or pheatmap will calculate the Z-scores from normalized counts/FPKM values? Could you please help me?
FPKM expression levels are not ideal for differential expression studies due to the fact that there is no cross-sample normalisation employed by the method that produces FPKM data. If you have already derived a list of statistically significantly differentially expressed genes, then I can only assume that you have done this via the results() function of DESeq2 (?), in which case, the test statistics can be thought to have been derived from the normalised counts.
In any case, you have 2 options:
transform your DESeq2 normalised counts via variance stabilisation
or regularised log (setting blind = FALSE, in either case), and
then directly running pheatmap on the transformed expression
levels, setting scale = 'row', i.e., pheatmap(..., scale =
'row'). This will scale the data to Z-scores (by gene)
automatically.
transform your FPKM expression levels to Z-scale directly via zFPKM
package, and then run pheatmap with scale = 'none' (as you will have
already scaled them to Z-scale via zFPKM). I show an example of this
with the old gplots::heatmap.2() on Biostars:
https://www.biostars.org/p/284721/#284723
I assume the rlog transformation has been done already(?) and I can use the normalized counts from the exported resdata directly for pheatmap if I set the scale = 'row' in pheatmap. Correct?
You create an object that stores regularised log expression levels with this line of code:
rld <- rlogTransformation(dds)
However, you never save that object or output it to disk (in the code snippet that you are showing)
I assume the rlog transformation has been done already(?) and I can
use the normalized counts from the exported resdata directly for
pheatmap if I set the scale = 'row' in pheatmap. Correct?
No, you should use assay(rld) as input to pheatmap(..., scale = 'row'). The dds object is independent from rld.
This isn’t a question about DESeq2 software per se, so I've removed the tag.