####material from //www2.warwick.ac.uk/fac/sci/moac/people/students/peter_cock/r/heatmap/ source("//www.bioconductor.org/biocLite.R") biocLite("ALL") library("ALL") data("ALL") ALL ##We can looks at the results of molecular biology testing for the 128 samples: ALL$mol.biol ### Ignoring the samples which came back negative on this test (NEG), ### most have been classified as having a translocation between chromosomes 9 and 22 (BCR/ABL), ### or a translocation between chromosomes 4 and 11 (ALL1/AF4). #### For the purposes of this example, we are only interested in these two subgroups, #### so we will create a filtered version of the dataset using this as a selection criteria: eset <- ALL[, ALL$mol.biol %in% c("BCR/ABL", "ALL1/AF4")] #### The resulting variable, eset, contains just 47 samples #### each with the full 12,625 gene expression levels. #### This is far too much data to draw a heatmap with, #### but we can do one for the first 100 genes as follows: heatmap(exprs(eset[1:100,])) ###use limma package to find differentially expressed genes source("//www.bioconductor.org/biocLite.R") biocLite("limma") library(limma) f <- factor(as.character(eset$mol.biol)) design <- model.matrix(~f) fit <- eBayes(lmFit(eset,design)) ## The leftmost numbers are row indices, ID is the Affymetrix HGU95av2 accession number, ## M is the log ratio of expression, A is the log average expression, t the moderated t-statistic, ## and B is the log odds of differential expression. topTable(fit, coef=2) ## Next, we select those genes that have adjusted p-values below 0.05, ## using a very stringent Holm method to select a small number (165) of genes. selected <- p.adjust(fit$p.value[, 2]) <0.05 esetSel <- eset [selected, ] ### now produce the heatmap heatmap(exprs(esetSel)) ### change the color scheme heatmap(exprs(esetSel), col=topo.colors(100)) ###we can use the color bars color.map <- function(mol.biol) { if (mol.biol=="ALL1/AF4") "#FF0000" else "#0000FF" } patientcolors <- unlist(lapply(esetSel$mol.bio, color.map)) heatmap(exprs(esetSel), col=topo.colors(100), ColSideColors=patientcolors) ##### One subtle point in the previous examples is that ##### the heatmap function has automatically scaled the colours for each row ##### (i.e. each gene has been individually normalised across patients). ##### scale = "none" heatmap(exprs(esetSel), col=topo.colors(75), scale="none", ColSideColors=patientcolors, cexRow=0.5) ####want to add color keys? ####use heatmap.2 in "gplots" library("gplots") heatmap.2(exprs(esetSel), col=topo.colors(75), scale="none", ColSideColors=patientcolors, key=TRUE, symkey=FALSE, density.info="none", trace="none", cexRow=0.5) ###green to red heatmap.2(exprs(esetSel), col=redgreen(75), scale="row", ColSideColors=patientcolors, key=TRUE, symkey=FALSE, density.info="none", trace="none", cexRow=0.5)