###################################################
### chunk number 1: setup
###################################################
options(width=60)
options(continue=" ")
options(prompt="R> ")
options(width=45, digits=2)


###################################################
### chunk number 2: loading
###################################################
  library(oligo)
  library(maqcExpression4plex)
  library(genefilter)
  library(limma)
  library(RColorBrewer)
  palette(brewer.pal(8, "Dark2"))


###################################################
### chunk number 3: listing
###################################################
  extdata <- system.file("extdata",
      package="maqcExpression4plex")
  xys.files <- list.xysfiles(extdata,
      full.names=TRUE)
  basename(xys.files)


###################################################
### chunk number 4: ExpressionFeatureSet
###################################################
  pd <- dir(extdata, pattern="phenoData", full.names=TRUE)
  pd <- read.AnnotatedDataFrame(pd)
  maqc <- read.xysfiles(xys.files, phenoData=pd)
  class(maqc)


###################################################
### chunk number 5: rawdata
###################################################
   exprs(maqc)[10001:10010, 1:2]


###################################################
### chunk number 6: maqcBoxplot
###################################################
  boxplot(maqc, main="MAQC Sample Data")


###################################################
### chunk number 7: maqcHist
###################################################
  hist(maqc, main="MAQC Sample Data")


###################################################
### chunk number 8: rma
###################################################
  eset <- rma(maqc)
  class(eset)
  show(eset)
  exprs(eset)[1:10, 1:2]


###################################################
### chunk number 9: rmaResultsB
###################################################
  boxplot(eset, transfo=identity, main="After RMA")


###################################################
### chunk number 10: rmaResultsH
###################################################
  hist(eset, transfo=identity, main="After RMA")


###################################################
### chunk number 11: naive
###################################################
  e <- exprs(eset)
  index <- which(eset[["Key"]] == "brain")
  d <- rowMeans(e[, index])-rowMeans(e[, -index])
  a <- rowMeans(e)
  sum(abs(d)>1)


###################################################
### chunk number 12: ttests
###################################################
  tt <- rowttests(e, factor(eset[["Key"]]))
  lod <- -log10(tt[["p.value"]])


###################################################
### chunk number 13: maplot
###################################################
  smoothScatter(a, d, xlab="Average Intensity", ylab="Log-ratio", main="MAQC Sample Data")
  abline(h=c(-1, 1), col=2)


###################################################
### chunk number 14: volcanoplot2
###################################################
  o1 <- order(abs(d),decreasing=TRUE)[1:25]
  o2 <- order(abs(tt[["statistic"]]),decreasing=TRUE)[1:25]
  o <- union(o1,o2)
  smoothScatter(d, lod, main="A Better view")
  points(d[o1], lod[o1], pch=18, col="blue")
  points(d[o2], lod[o2], pch=1, col="red")
  abline(h=2, v=c(-1, 1), col=2)


###################################################
### chunk number 15: limma
###################################################
  design <- model.matrix(~factor(eset[["Key"]]))
  fit <- lmFit(eset, design)
  ebayes <- eBayes(fit)
  lod <- -log10(ebayes[["p.value"]][,2])
  mtstat<- ebayes[["t"]][,2]


###################################################
### chunk number 16: volcanoplot3
###################################################
  o1 <- order(abs(d), decreasing=TRUE)[1:25]
  o2 <- order(abs(mtstat), decreasing=TRUE)[1:25]
  o <- union(o1, o2)
  smoothScatter(d, lod, main="Moderated t", xlab="Log-ratio", ylab="LOD")
  points(d[o1], lod[o1], pch=18,col="blue")
  points(d[o2], lod[o2], pch=1,col="red")
  abline(h=2, v=c(-1, 1))


###################################################
### chunk number 17: toptable
###################################################
  tab <- topTable(ebayes, coef=2, adjust="fdr", n=10)
  tab


###################################################
### chunk number 18: sessionInfo
###################################################
  sessionInfo()