###################################################
### chunk number 1: 
###################################################
library(affycoretools)
library(KEGG)
library(xtable)
library(affyPLM)
## make AnnotatedDataFrame
pd <- read.AnnotatedDataFrame(paste(system.file("examples",package="affycoretools"),
                                    "/pdata.txt", sep=""), header = TRUE, row.names = 1)

## no celfiles in package any more, fake this step
#dat <- ReadAffy()
#eset <- rma(dat)

load(paste(system.file("examples", package="affycoretools"),
           "/abatch.Rdata", sep=""))
load(paste(system.file("examples", package="affycoretools"),
           "/exprSet.Rdata", sep=""))

## load annotation package
options(show.error.messages = FALSE)
a <- try(do.call("library", list(annotation(eset))))
options(show.error.messages = TRUE)
if(inherits(a, "try-error")){
source("http://www.bioconductor.org/biocLite.R")
biocLite(annotation(eset))
do.call("library", list(annotation(eset)))
}
type <- rep(LETTERS[1:4], each = 3)


###################################################
### chunk number 2: 
###################################################
## Put a table of filenames and sample types in PDF
type <- rep(LETTERS[1:4], each = 3)
tmp <- data.frame(sampleNames(eset), type)
names(tmp) <- c("Filenames","Samples")
xtable(tmp)


###################################################
### chunk number 3: 
###################################################
plotHist(dat, sampleNames(eset))
plotHist(dat[,1:6])
plotHist(dat[,7:12])


###################################################
### chunk number 4: 
###################################################
plotDeg(dat, sampleNames(eset))
plotDeg(dat[,1:6])
plotDeg(dat[,7:12])


###################################################
### chunk number 5: 
###################################################
boxplot(data.frame(exprs(eset)))


###################################################
### chunk number 6: 
###################################################
plotPCA(eset, groups = rep(1:4, each = 3), 
        groupnames = unique(paste(pData(pd)[,1], pData(pd)[,2], sep = "-")))
plotPCA(eset[,1:6], groups=rep(1:2, each=3),
        groupnames=unique(paste(pData(pd)[1:4,1], pData(pd)[1:4,2], sep="-")))
plotPCA(eset[,7:12], groups=rep(1:2, each=3),
        groupnames=unique(paste(pData(pd)[7:10,1], pData(pd)[7:10,2], sep="-")))



###################################################
### chunk number 7: 
###################################################
library(affyPLM)
pset <- fitPLM(dat[,7:12])
boxplot(pset, xaxt="n", main="NUSE")
axis(1, at=7:12, label=paste("Sample", 1:6), las=2, cex.axis=0.8)


###################################################
### chunk number 8: 
###################################################
Mbox(pset, xaxt="n", main="RLE")
axis(1, at=1:6, label=paste("Sample", 1:6), las=2, cex.axis=0.8)


###################################################
### chunk number 9: 
###################################################
## filter data to remove probesets that aren't changing
index <- apply(exprs(eset)[,1:6], 1, var) > 0.01
eset1 <- eset[index,]
## create design matrix and give reasonable column names
design <- model.matrix(~ 0 + factor(rep(1:4,each=3)))
colnames(design) <- LETTERS[1:4]


###################################################
### chunk number 10: 
###################################################
## set up a contrasts matrix using makeContrasts()
contrast <- makeContrasts(A - B, C - D, levels = design)

## now do the same using matrix()
contrast <- matrix(c(1,-1,0,0,0,0,1,-1), ncol=2, 
                   dimnames=list(unique(type), paste(type[c(1,7)],type[c(4,10)],
                   sep=" vs ")))


###################################################
### chunk number 11: 
###################################################
## fit model
fit <- lmFit(eset1, design)
fit2 <- contrasts.fit(fit, contrast)
## empirical Bayes step
fit2 <- eBayes(fit2)


###################################################
### chunk number 12: 
###################################################
## output text and HTML tables using limma2annaffy()
out <- limma2annaffy(eset1, fit2, design, contrast, annotation(eset),
                     pfilt = 0.05, fldfilt = 1, save = TRUE, text = TRUE,
                     interactive = FALSE)


###################################################
### chunk number 13: 
###################################################
## output text and HTML tables using limma2biomaRt()
## this is usually more applicable to analyses of chips that don't have BioC 
## annotation packages
## Change mysql to FALSE if you don't have RMySQL installed

out2 <- limma2biomaRt(eset1, fit2, design, contrast, "hsapiens",
                      ann.source="affy_hg_focus", pfilt=0.05, fldfilt=1, 
                      addname="biomart", affyid=TRUE, mysql=TRUE, interactive=FALSE)


###################################################
### chunk number 14: 
###################################################
a <- data.frame(paste(unique(type)[c(1,3)], unique(type)[c(2,4)], sep = " vs "),
                sapply(out, function(x) dim(x)[1]))
names(a) <- c("Comparisons","Probesets")
xtable(a, digits = rep(0, 3))