Select a list of genes of interest to visualise in heatmap of DESeq object

Question

[jason.williams]

Hi,

I have previously generated heatmaps of the top x genes in DESeq2 using the following code:

topVarGenes <- head(order(rowVars(assay(rld)), decreasing = TRUE), 20)
mat  <- assay(rld)[ topVarGenes, ]
mat  <- mat - rowMeans(mat)
anno <- as.data.frame(colData(rld)[c("Condition")])
pheatmap(mat, annotation_col = anno)

This is in the main DESeq2 vignette and a heat map is generated with conditions and patients metadata etc.

As this just selects the top x genes, I was wondering if I can also select a list of genes of interest by ensembl ID, so that I can generate a heat map to show how these genes of interest change between conditions. This would also be useful to see how a list of genes associated with an ontology change between conditions.

Best wishes,

Jason

Answer 1

[rpolicastro]

If ensembl IDs were used to generate the count matrix, you can simply replace topVarGenes with a vector that contains any genes you wish to test.

genes <- c("gene_1", "gene_2", ...)
mat <- assay(rld)[genes, ]
mat <- mat - rowMeans(mat)
anno <- as.data.frame(colData(rld)[c("Condition")])
pheatmap(mat, annotation_col = anno)

[cpad0112]

mat <- mat - rowMeans(mat) what is the function of this line? @ rpolicastro

[rpolicastro]

This was a line in their original code that I carried over. For each gene, the expression value for that gene in each sample is subtracted by the mean expression value for that gene over all samples.It's a way to center the data for display in a heatmap.

[cpad0112]

pheatmap has in-built row-wise or column-wise scaling option. One can use that function assuming that user wants to center the data row wise. On the other hand, i think z-scaling is more useful than subtracting means from the values (subtract and divide by standard deviation).

[jason.williams]

This worked really nicely, thanks very much @rpolicastro, I think this will be useful to others that want to select a custom set of genes of interest in an RNA-Seq experiment.

[abrilluchetti]

Hi! I'm trying to make a pheatmap with a group of 600 genes. When I make the list like you did, I get an error because of the lenght of the code. I tried dividing the list into 5 vectors and then using c() to make one, but I get the error "subscript out of bounds". When I import the list from an excel file and use it as a data frame I get this error: invalid subscript type 'list'. I have tried to change this list to different types but I don't understand what's wrong with it. Do you have any idea of what could be happening? Thanks!!

[Pankoj]

What is assay here? It gives error "could not find function "assay"

[rpolicastro]

assay gets the count matrix from SummarizedExperiment or related objects like the DESeq2 object.