###################################################
### chunk number 1: setup
###################################################
  options(width=45, digits=2)


###################################################
### chunk number 2: loading
###################################################
  library(oligo)
  library(maqcExpression4plex)
  library(genefilter)
  library(geneplotter)
  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
###################################################
  maqc <- read.xysfiles(xys.files)
  pd <- dir(extdata, pattern="phenoData", full.names=TRUE)
  phenoData(maqc) <- read.AnnotatedDataFrame(pd)
  class(maqc)


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


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


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


###################################################
### chunk number 10: rmaResults
###################################################
  hist(eset, main="After RMA")


###################################################
### chunk number 11: fix
###################################################
  options(width=35)


###################################################
### chunk number 12: naive
###################################################
  e <- exprs(eset)
  dim(e)
  index <- 1:3
  d <- rowMeans(e[, index])-rowMeans(e[, -index])
  a <- rowMeans(e)
  sum(abs(d)>1)


###################################################
### chunk number 13:  eval=FALSE
###################################################
##   smoothScatter(a, d, xlab="Average Intensity", ylab="Log-ratio", main="MAQC Sample Data")


###################################################
### chunk number 14: naiveMAplot
###################################################
  smoothScatter(a, d, xlab="Average Intensity", ylab="Log-ratio", main="MAQC Sample Data")


###################################################
### chunk number 15: ttests
###################################################
  tt <- rowttests(e, factor(eset$Key))
  lod <- -log10(tt$p.value)


###################################################
### chunk number 16:  eval=FALSE
###################################################
##   smoothScatter(d, lod, xlab="Log-ratio", ylab="LOD", main="MAQC Sample Data")
##   abline(h=2, v=c(-1, 1))


###################################################
### chunk number 17: ttestsPlot
###################################################
  smoothScatter(d, lod, xlab="Log-ratio", ylab="LOD", main="MAQC Sample Data")
  abline(h=2, v=c(-1, 1))


###################################################
### chunk number 18: 
###################################################
  o1 <- order(abs(d),decreasing=TRUE)[1:25]
  o2 <- order(abs(tt$statistic),decreasing=TRUE)[1:25] #$
  o <- union(o1,o2)


###################################################
### chunk number 19:  eval=FALSE
###################################################
##   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")


###################################################
### chunk number 20: volcanoplot2
###################################################
  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")


###################################################
### chunk number 21: fix
###################################################
  options(width=45)


###################################################
### chunk number 22: 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 23: 
###################################################
  o1 <- order(abs(d),decreasing=TRUE)[1:25]
  o2 <- order(abs(mtstat),decreasing=TRUE)[1:25]
  o <- union(o1,o2)


###################################################
### chunk number 24:  eval=FALSE
###################################################
##   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 25: volcanoplot3
###################################################
  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 26: fix
###################################################
  options(width=60)


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


###################################################
### chunk number 28: fix
###################################################
  options(width=35)


###################################################
### chunk number 29: prepare
###################################################
library("oligo")
library("hapmap100kxba")
pathCelFiles <- system.file("celFiles", package="hapmap100kxba")
fullFilenames <- list.celfiles(path=pathCelFiles, full.names=TRUE)


###################################################
### chunk number 30: getSQS
###################################################
temporaryDir <- tempdir()
preProcessedData <- justSNPRMA(fullFilenames, tmpdir=temporaryDir)
preProcessedData$gender <- c("female", "female", "male")


###################################################
### chunk number 31: 
###################################################
theA <- getA(preProcessedData)
theM <- getM(preProcessedData)
dim(theA)


###################################################
### chunk number 32:  eval=FALSE
###################################################
## smoothScatter(theA[,1,1],
##               theM[,1,1],
##               main="MA-plot (Antisense)",
##               xlab="Average Intensity",
##               ylab="Log-ratio (A/B)")


###################################################
### chunk number 33: SNPmaplot
###################################################
smoothScatter(theA[,1,1], theM[,1,1], main="MA-plot (Antisense)",
              xlab="Average Intensity", ylab="Log-ratio (A/B)")


###################################################
### chunk number 34: crlmm
###################################################
crlmmOut <- crlmm(preProcessedData,
                  correctionFile="exampleCorrection.rda",
                  verbose=FALSE)
calls(crlmmOut)[1:3,1]
range(callsConfidence(crlmmOut))


###################################################
### chunk number 35: fix
###################################################
  options(width=40)


###################################################
### chunk number 36: sql1
###################################################
conn <- db(crlmmOut)
dbListTables(conn)


###################################################
### chunk number 37: sql2
###################################################
dbListFields(conn, "featureSet")


###################################################
### chunk number 38: sql3
###################################################
fields <- c("man_fsetid, chrom, physical_pos")
cond <- c("man_fsetid LIKE 'SNP%' LIMIT 5")
sql <- paste("SELECT", fields, "FROM featureSet WHERE", cond)
dbGetQuery(conn, sql)


###################################################
### chunk number 39: sql4
###################################################
p1 <- "SELECT man_fsetid, physical_pos"
p2 <- "FROM featureSet WHERE man_fsetid"
p3 <- "LIKE 'SNP%' AND chrom='X'"
p4 <- "ORDER BY physical_pos"
sql <- paste(p1, p2, p3, p4)
x.info <- dbGetQuery(conn, sql)


###################################################
### chunk number 40: 
###################################################
  idx <- match(x.info[,1], rownames(theA))
tmpA <- rowMeans(theA[idx,,], dims=2)


###################################################
### chunk number 41:  eval=FALSE
###################################################
## plot(1, type="n", xlab="Physical Position", ylab="Average Intensity",
##      main="Intensities on Chromosome X", ylim=c(10.5, 12),
##      xlim=range(x.info[, 2]))
## for (i in 1:3)
##   lines(lowess(x.info[, 2], tmpA[, i]), col=i, lwd=2)
## legend("top", paste("Sample ", 1:3), col=1:3, lwd=2, lty=1)


###################################################
### chunk number 42: xchr2
###################################################
plot(1, type="n", xlab="Physical Position", ylab="Average Intensity",
     main="Intensities on Chromosome X", ylim=c(10.5, 12),
     xlim=range(x.info[, 2]))
for (i in 1:3)
  lines(lowess(x.info[, 2], tmpA[,i]), col=i, lwd=2)
legend("top", paste("Sample ", 1:3), col=1:3, lwd=2, lty=1)


###################################################
### chunk number 43: 
###################################################
  options(width=60)


###################################################
### chunk number 44: sessionInfo
###################################################
  sessionInfo()