Heatmap for gene cluster

Question

[dazhudou1122]

Dear Biostarts community,

I am drawing a heatmap for gene clustering based on row variation, using the results from a standard DEseq2 pipeline.

I used the following code to pick top varied genes:

topVarGenes <- head(order(rowVars(assay(rld)), decreasing = TRUE), 100)

but the problems is that the head function skewed the results. I then changed the code to

topVarGenes <- rld[order(rowVars(assay(rld))), decreasing = TRUE, ]

then I got error:

Error in slot(object, slot_name) :
  no slot of name "decreasing" for this object of class "DESeqTransform"

Can someone help me figure out what went wrong? I would very much appreciate that.

Thank you!

Wenhan

I have put the complete code below:

# pheatmap for gene clustering


library("genefilter")
  library(pheatmap)
  topVarGenes <- head(order(rowVars(assay(rld)), decreasing = TRUE), 100)
  rld <- rld[order(rowVars(assay(rld)),decreasing = TRUE),]
  topVarGenes <- rld[1:100,]
  mat  <- assay(rld)[ topVarGenes, ]
  mat  <- mat - rowMeans(mat)
  anno <- as.data.frame(colData(rld))
  row.names(mat) <- gsub('\\..+$', '', row.names(mat))
  library(org.Hs.eg.db)
  gns <- select(org.Hs.eg.db, row.names(mat), "SYMBOL", "ENSEMBL")
  row.names(mat)[match(gns[,1], row.names(mat))] <- gns[,2]
  png("gene_clustering_pheatmap.png", 1500, 3000, pointsize=25)
  pheatmap(mat, annotation_col = anno, fontsize=24)
  dev.off()

Answer 1

[Kevin Blighe]

Hi, you just need to do this, for selecting the top 100 most variable genes:

library(airway)
library(magrittr)

data('airway')
airway$dex %<>% relevel('untrt')

library('DESeq2')

dds <- DESeqDataSet(airway, design = ~ cell + dex)
dds <- DESeq(dds, betaPrior=FALSE)
rld <- rlog(dds)

top100 <- head(order(rowVars(assay(rld)), decreasing = TRUE), 100)
mat <- as.matrix(assay(rld)[top100,])

[dazhudou1122]

Thank you Kevin for your help. This is a very clean and elegant solution. It is working properly now.