####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)