Golub dataset analysis:

calc.p.value <- function(v, x, y) {  
  if ( sd(v[c(x,y)]) == 0 ) {
    return(1)
  } else {
    r <- t.test(v[x], v[y])
    return(r$p.value)
  }
}

golub <- read.delim("http://www.bioinf.ucd.ie/people/ian/Golub.txt", sep="\t")
head(golub)
colnames(golub)
library(limma)
samples.AML <- colnames(golub)[ grep("^AML", colnames(golub)) ]
samples.ALL <- colnames(golub)[ grep("^ALL", colnames(golub)) ]
golub.samps <- c(samples.AML, samples.ALL)
golub[, golub.samps] <- normalizeQuantiles(golub[,c(samples.AML,samples.ALL)])


golub$ALL.mean <- apply(golub[, samples.ALL], 1, mean)
golub$AML.mean <- apply(golub[, samples.AML], 1, mean)


golub$logFC <- log2(golub$AML.mean/golub$ALL.mean)


golub$p.value <- apply(golub[,c(samples.AML,samples.ALL)], 1, calc.p.value, samples.AML, samples.ALL)
golub$q.value <- p.adjust( golub$p.value , method="fdr" )



# AML up-regulated
dim( subset(golub, q.value<0.05 & logFC > 0.5) )

# ALL up-regulated
dim( subset(golub, q.value<0.05 & logFC < -0.5) )

# TOGETHER
dim( subset(golub, q.value<0.05 & (logFC > 0.5 | logFC < -0.5)) )

dim( subset(golub, q.value<0.05 & abs(logFC) > 0.5 ) )

df.golub.diff <- subset(golub, q.value<0.05 & abs(logFC) > 0.5 )


golub.association <- read.delim("/tmp/golub_association.csv",sep=",")

colnames(golub.association)


df.golub.diff.annot <- merge(x=df.golub.diff,
      y=golub.association,
      by.x='Golubtemplate',
      by.y='Gene.Accession.Number'
      )

write.table(df.golub.diff[,c('Golubtemplate','AML.mean','ALL.mean','logFC','p.value','q.value')], file="/tmp/gobub_diff.tsv", 
            col.names=TRUE, row.names=FALSE, quote=FALSE, sep="\t"
            )
boxplot(df.golub.diff[,c('AML.mean','ALL.mean')])



library(gplots)


#svg(filename="/tmp/golub.diffs.svg", 
#    width = 8, height = 6)


png(filename = "/tmp/golub.diffs.png", 
    bg="white", 
    res=300, 
    width=3000, height=3000)

hv1<-heatmap.2(
  as.matrix(log2(df.golub.diff.annot[,golub.samps])),
  scale="row",
  col=greenred,
  labRow=df.golub.diff.annot[,"Gene.Description"],
  key=TRUE,
  symkey=TRUE,
  density.info="none",
  trace="none",
  Rowv=T,
  Colv=T,
  cexRow=0.45,
  cexCol=1,
  keysize=1,
  margins = c(10,25),
  dendrogram=c("both"),
  main = "HeatMap"
)

graphics.off()