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
•
written
by
megannj •
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.
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.
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).
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
Log in to answer this question.
