Question: Plot LFC with pheatmap of differentially expressed gene list from DESeq2.
1
gravatar for megannj
9 months ago by
megannj10
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
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

deseq deseq2 pheatmap heatmap • 834 views
ADD COMMENTlink modified 9 months ago by Devon Ryan89k • written 9 months ago by megannj10
1
gravatar for Devon Ryan
9 months ago by
Devon Ryan89k
Freiburg, Germany
Devon Ryan89k 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 9 months ago • written 9 months ago by Devon Ryan89k

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 REPLYlink written 9 months ago by megannj10
1

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 9 months ago by Devon Ryan89k
Please log in to add an answer.

Help
Access

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