## ----include=FALSE-------------------------------------------------------
knitr::opts_chunk$set(tidy=FALSE)


## ----style, eval=TRUE, echo=FALSE, results="asis"------------------------
BiocStyle::latex()


## ----eval = FALSE--------------------------------------------------------
## library("affycoretools")
## pd <- new("AnnotatedDataFrame",
##           data = read.table("pdata.txt", header = TRUE, row.names = 1))
## eset <- affystart(groups = rep(1:4, each = 3),
##                   groupnames = unique(paste(pData(pd)[,1],
##                   pData(pd)[,2], sep = "-")),
##                   phenoData = pd)


## ----echo = FALSE, include=FALSE-----------------------------------------
library(affycoretools)
load("abatch.Rdata")
load("exprSet.Rdata")


## ----hist, fig.cap = "Density plot", echo = FALSE, fig.width = 5, fig.height = 5----
plotHist(dat)


## ----deg, fig.cap = "RNA degradation plot", echo = FALSE, fig.width = 5, fig.height = 5----
plotDeg(dat)


## ----pca, fig.cap = "PCA plot", echo = FALSE, fig.width = 5, fig.height = 5----
pd <- new("AnnotatedDataFrame", 
          data = read.table("pdata.txt", header = TRUE, row.names = 1))
sampleNames(pd) <- sampleNames(eset)
plotPCA(eset, groups = rep(1:4, each = 3), 
        groupnames = unique(paste(pData(pd)[,1], pData(pd)[,2], sep = "-")))
phenoData(eset) <- pd


## ----scree, fig.cap = "Screeplot", echo = FALSE, fig.width = 5, fig.height = 5----
plotPCA(eset, screeplot = TRUE)


## ----eval = FALSE--------------------------------------------------------
## eset1 <- affystart(filenames = list.celfiles()[1:6],
##                    plot = FALSE, pca = FALSE)
## eset2 <- affystart(filenames = list.celfiles()[7:12],
##                    plot = FALSE, pca = FALSE)
## eset <- new("ExpressionSet",
##             exprs = cbind(exprs(eset1), exprs(eset2)),
##             phenoData = pd,
##             annotation = annotation(eset1))


## ------------------------------------------------------------------------
library(genefilter)
f1 <- kOverA(3, 6)
filt <- filterfun(f1)
index <- genefilter(eset, filt)
eset <- eset[index,]


## ------------------------------------------------------------------------
library(limma)
grps <- paste(pData(eset)[,1], 
              pData(eset)[,2], sep = ".")
design <- model.matrix(~ 0 + factor(grps))
colnames(design) <- levels(factor(grps))
ugrps <- unique(grps)
contrasts <- matrix(c(1, -1, 0, 0, 0, 0, 1, -1),
                    ncol = 2, dimnames = list(ugrps,
                              paste(ugrps[c(1,3)], ugrps[c(2,4)],
                                    sep = " - ")))
fit <- lmFit(eset, design)
fit2 <- contrasts.fit(fit, contrasts)
fit2 <- eBayes(fit2)


## ------------------------------------------------------------------------
design
contrasts


## ----venn, fig.cap = "Venn Diagram"--------------------------------------
rslt <- decideTests(fit2)
vc <- vennCounts2(rslt, method = "same")
vennDiagram(vc, cex = 0.8)


## ----eval = FALSE--------------------------------------------------------
## vennSelect(eset, design, rslt, contrasts, fit2)


## ----eval = FALSE--------------------------------------------------------
## limma2annaffy(eset, fit2, design,
##               contrasts, annotation(eset),
##               pfilt = 0.05)


## ----eval = FALSE--------------------------------------------------------
## index1 <- vennSelect(x = rslt, indices.only = TRUE)[[3]]
## probids <- unique(getLL(featureNames(eset)[index1],
##                         annotation(eset)))
## univ <- unique(getLL(featureNames(eset),
##                      annotation(eset)))
## params <- new("GOHyperGParams", geneIds = probids,
##               universeGeneIds = univ,
##               annotation = annotation(eset),
##               conditional = TRUE, ontology = "MF")
## hyp <- hyperGTest(params)
## htmlReport(hyp, file = "GO MF terms.html",
##            categorySize = 10)


## ----sessioninfo, echo=FALSE, results="markup", comment=NA---------------
sessionInfo()