Plot LFC with pheatmap of differentially expressed gene list from DESeq2.

Question

[megannj]

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
source("https://bioconductor.org/biocLite.R")
biocLite("DESeq2")
biocLite("stringi")
biocLite("MASS")
install.packages("survival")
biocLite("bit")
library("DESeq2")
setwd("C:/Users/megan/Desktop/Computational/WizKO-RNAseq")

install.packages("rlang")
biocLite("DESeqDataFromMatrix") 
##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)
resultsNames(dds)

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

   ##Heatmap 
biocLite("pheatmap")
library("pheatmap")
ntd <- normTransform(dds)
select <- order(rowMeans(counts(dds, normalized=TRUE)), decreasing =TRUE)[1:50]
df <- as.data.frame(colData(dds)[,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.

Megan

Answer 1

[Devon Ryan]

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.

[megannj]

Devon,

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.

[Devon Ryan]

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).

[margott.j]

Hi Megan, this was some time ago but I was wondering if you ever figured out how to manipulate the code to plot the LFC? I have tried to figure this out for a while now and have been unsuccessful so far.

Thank you!

Margot