Question: Plot LFC with pheatmap of differentially expressed gene list from DESeq2.
gravatar for megannj
2.0 years ago by
megannj10 wrote:

Hi, all!

First post, so apologies for any flaws with post structure.

I am attempting to make a basic heatmap that shows the log fold change of differentially expressed genes, as identified by DESeq2. See below the code I am using for DESeq2:

##Load DESeq2

##Read in your data
Counts_bbmap <- read.table("RNAseqResults.txt", header = TRUE, row.names=1)
sampleInfo <- read.csv("ValuesFile.txt")
sampleInfo <- data.frame(sampleInfo)
dds <- DESeqDataSetFromMatrix(countData = Counts_bbmap, colData = sampleInfo, design = ~Genotype)
dds <- DESeq(dds)

To make my heatmap, I am attempting to use pheatmap using the following code:

ntd <- normTransform(dds)
select <- order(rowMeans(counts(dds, normalized=TRUE)), decreasing =TRUE)[1:50]
df <-[,c("Genotype","Sample")])
pheatmap(assay(ntd)[select,], cluster_rows=FALSE, showrownames=TRUE, cluster_cols=FALSE, annotation_col=df)

However, I want to only post Log Fold Changes, and the above code produces a heatmap plotting something else (although I'm not sure what). I've tried a few different things myself, such as subsetting the data into just genes with adjusted p values over a certain number, sorting the data, etc. But I cannot figure out how to just plot the log fold changes that have already been calculated. The error may be in the assay(ntd) but I do not know how to resolve this.

Thanks for any and all help.


deseq deseq2 pheatmap heatmap • 2.2k views
ADD COMMENTlink modified 2.0 years ago by Devon Ryan95k • written 2.0 years ago by megannj10
gravatar for Devon Ryan
2.0 years ago by
Devon Ryan95k
Freiburg, Germany
Devon Ryan95k wrote:

Instead of plotting assay(ntd), you'll want to instead plot lfcShrink(dds)$log2FoldChange. Realistically, use something like:

res = lfcShrink(dds)
... whatever filtering you want can go here ...
pheatmap(res$log2FoldChange, ...)

Now having written that, I would suggest that you NOT make a heatmap of fold-changes. I have yet to see a case where these were ever actually informative. I think the heatmap in your original code is MUCH more useful, since you can actually see how genes might relate to each other in terms of covariation.

ADD COMMENTlink modified 2.0 years ago • written 2.0 years ago by Devon Ryan95k


Thank you so much! I will try this! I think the heatmap I've made with the code above would be good if I could filter or sort the output somehow. Right now it is cluttered by genes that are not expressed in our cell type, so a lot of the rows are just solid blue all the way across.

ADD REPLYlink written 2.0 years ago by megannj10

Instead of plotting the top 50 genes by mean signal (what you're doing now), you might try plotting be top variance instead. I think there's an example of this in the DESeq2 vignette (if not, there are a number of nice QC plots in there that you can use for inspiration).

ADD REPLYlink written 2.0 years ago by Devon Ryan95k
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 807 users visited in the last hour