Plot LFC with pheatmap of differentially expressed gene list from DESeq2.
1
2
Entering edit mode
6.3 years ago
megannj ▴ 20

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

DESeq2 DESeq pheatmap heatmap • 5.9k views
ADD COMMENT
1
Entering edit mode
6.3 years ago

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 COMMENT
0
Entering edit mode

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.

ADD REPLY
1
Entering edit mode

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 REPLY
0
Entering edit mode

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

ADD REPLY

Login before adding your answer.

Traffic: 1949 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6