How to make heatmap of log fold change for multiple conditions from RNA-seq

Question

[mar77]

I'm trying to make a heatmap of the most significantly differentially expressed genes from some RNA-seq data. I have several different conditions and a control condition with triplicate repeats for each condition. I want to plot the log fold change for the different samples to visualise how they differ from the control rather than how they differ from the average across all samples (which is what I believe heatmaps normally show and I have been able to generate) - so as the control is no change it wouldn't necessarily be included in the plot.

So far I have done:

ddsObj.raw<- DESeqDataSetFromMatrix(countData=data, colData=sampleinfo, design=~Condition)
ddsObj<-DESeq(ddsObj.raw)

I know I can use

res<- results(ddsObj, contrast=c("Condition","Treatment1","Control")
log2_changes<-res[,"log2FoldChange"]

to extract the log2 fold changes for a specified pair wise comparison.

Is there a way I can extend this to include the other conditions as well, i.e. to also have treatment2 v control, treatment3 v control, etc and to use this to plot a heatmap of the most significantly differentially expressed genes. Is there a way to do this that is easier than going through all the pairwise comparisons individually and merging them together into one data frame?

I'm relatively new to bioinformatics so I hope this is ok to ask, if anyone can help or point me in the direction of a good tutorial it would be really, really appreciated!

Answer 1

[Kevin Blighe]

Hi, this section in the vignette explains it quite well: Contrasts.

Note that you should also be performing fold-change shrinkage, and, to avail of the contrast functionality with fold-change shrinkage, you should use type = 'ashr' (see Extended section on shrinkage estimators).

For specifying different comparisons via contrast, I gave some examples here: A: DESeq2 compare all levels

So, for all intents and purposes, you'd need [for treatment2 versus Control]:

res <- results(ddsObj, contrast = c('Condition', 'Treatment2', 'Control'))
res <- lfcShrink(dds, contrast = c('Condition', 'Treatment2', 'Control'),
  res = res, type = 'ashr')

########

For the heatmap part, you can use the subset() command to subset your results object for statistically significantly differently expressed genes. You then need to derive transformed expression levels from DESeq2, preferentially via the vst() function.

Once you have a list of statistically significantly differently expressed genes and a transformed expression matrix, you need subset this matrix for these genes and then generate the heatmap on this subsetted matrix, something like:

vsd <- vst(dds)
heat <- assay(vsd)[GeneList,]

I have a heatmap tutorial here that requires no further coding on your part, but it is advanced an would require adaptation for your own data: A simple tutorial for a complex ComplexHeatmap

Kevin