Heatmaps in R

Although not the most interactive tool, R is great for generating plots of data very quickly.  We still use heatmaps a lot here for gene expression data, and so for a bit of nostalgia, I have some code below that draws a heatmap of the classic Yeast cell cycle data by Spellman et al.  This was a seminal microarray paper from Stanford!  This should execute in R:

# You will need to install the gplots and geneplotter libraries:


# Mol Biol Cell. 1998 Dec;9(12):3273-97.
# Comprehensive identification of cell cycle-regulated genes of the yeast
# Saccharomyces cerevisiae by microarray hybridization.
# Spellman PT, Sherlock G, Zhang MQ, Iyer VR, Anders K, Eisen MB, Brown PO, Botstein D, Futcher B.
# Department of Genetics, Stanford University Medical Center, Stanford, California 94306-5120, USA.
# We sought to create a comprehensive catalog of yeast genes whose transcript levels vary
# periodically within the cell cycle. To this end, we used DNA microarrays and samples from
# yeast cultures synchronized by three independent methods: alpha factor arrest, elutriation,
# and arrest of a cdc15 temperature-sensitive mutant. Using periodicity and correlation algorithms,
# we identified 800 genes that meet an objective minimum criterion for cell cycle regulation.
# In separate experiments, designed to examine the effects of inducing either the G1 cyclin
# Cln3p or the B-type cyclin Clb2p, we found that the mRNA levels of more than half of these
# 800 genes respond to one or both of these cyclins. Furthermore, we analyzed our set of cell
# cycle-regulated genes for known and new promoter elements and show that several known elements
# (or variations thereof) contain information predictive of cell cycle regulation. A full
# description and complete data sets are available at http://cellcycle-www.stanford.edu

# read the data in from URL
bots <- read.table(url("http://genome-www.stanford.edu/cellcycle/data/rawdata/combined.txt"), sep="\t", header=TRUE)

# get just the alpha data
abot <- bots[,c(8:25)]
rownames(abot) <- bots[,1]

# get rid of NAs
abot[is.na(abot)] <- 0

# we need to find a way of reducing the data. Can't do ANOVA as there are no
# replicates. Sort on max difference and take first 1000

min <-apply(abot, 1, min)
max <- apply(abot, 1, max)
sabot <- abot[order(max - min, decreasing=TRUE),][1:1000,]

# cluster on correlation
cdist <- as.dist(1 - cor(t(sabot)))
hc <- hclust(cdist, "average")

# draw a heatmap



And this is the result: