### R code from vignette source 'vignettes/beadarray/inst/doc/beadsummary.rnw'

###################################################
### code chunk number 1: installExampleData (eval = FALSE)
###################################################
## source("http://www.bioconductor.org/biocLite.R")
## biocLite(c("beadarrayExampleData", "illuminaHumanv3.db"))
## 


###################################################
### code chunk number 2: prelim
###################################################
library("beadarray")

require(beadarrayExampleData)

data(exampleSummaryData)

exampleSummaryData




###################################################
### code chunk number 3: objectPeek
###################################################

exprs(exampleSummaryData)[1:5,1:5]
se.exprs(exampleSummaryData)[1:5,1:5]




###################################################
### code chunk number 4: annotations
###################################################

head(fData(exampleSummaryData))
table(fData(exampleSummaryData)[,"Status"])

pData(exampleSummaryData)



###################################################
### code chunk number 5: subsetting1
###################################################

channelNames(exampleSummaryData)

exampleSummaryData.log2 <- channel(exampleSummaryData, "G")
exampleSummaryData.unlogged <- channel(exampleSummaryData, "G.ul")


sampleNames(exampleSummaryData.log2)
sampleNames(exampleSummaryData.unlogged)

exprs(exampleSummaryData.log2)[1:10,1:3]
exprs(exampleSummaryData.unlogged)[1:10,1:3]



###################################################
### code chunk number 6: subsetting2
###################################################


exampleSummaryData.log2[,1:4]
exampleSummaryData.log2[1:10,]



###################################################
### code chunk number 7: subsetting4
###################################################

randIDs <- sample(featureNames(exampleSummaryData), 1000)

exampleSummaryData[randIDs,]



###################################################
### code chunk number 8: boxplot1
###################################################

boxplot(exampleSummaryData.log2[randIDs,])



###################################################
### code chunk number 9: boxplot2
###################################################

boxplot(exampleSummaryData.log2[randIDs,], what="nObservations")



###################################################
### code chunk number 10: boxplot4
###################################################

boxplot(exampleSummaryData.log2[randIDs,], sampleFactor="SampleFac")



###################################################
### code chunk number 11: boxplot5
###################################################

boxplot(exampleSummaryData.log2[randIDs,], probeFactor = "Status")



###################################################
### code chunk number 12: addFdata
###################################################

annotation(exampleSummaryData)

exampleSummaryData.log2 <- addFeatureData(exampleSummaryData.log2, toAdd = c("SYMBOL", "PROBEQUALITY", "CODINGZONE", "PROBESEQUENCE", "GENOMICLOCATION"))

head(fData(exampleSummaryData.log2))

illuminaHumanv3()



###################################################
### code chunk number 13: boxplot6
###################################################
ids <- which(fData(exampleSummaryData.log2)[,"SYMBOL"] == "ALB")

boxplot(exampleSummaryData.log2[ids,], sampleFactor = "SampleFac", probeFactor = "IlluminaID")


###################################################
### code chunk number 14: ggplot-layout
###################################################
require("gridExtra")
bp1 <- boxplot(exampleSummaryData.log2[ids,], sampleFactor = "SampleFac", probeFactor = "IlluminaID") + labs(title = "ALB expression level comparison") + xlab("Illumina Probe") + ylab("Log2 Intensity")

bp2 <- boxplot(exampleSummaryData.log2[randIDs,], probeFactor = "Status") + labs(title = "Control Probe Comparison")

print(bp1, vp = viewport(width = 0.5, height = 1, x = 0.25, y = 0.5))
print(bp2, vp = viewport(width = 0.5, height = 1, x = 0.75, y = 0.5))





###################################################
### code chunk number 15: MAs
###################################################
plotMAXY(exprs(exampleSummaryData), arrays=1:3,pch=16,log=FALSE)


###################################################
### code chunk number 16: normalise1
###################################################

exampleSummaryData.norm <- normaliseIllumina(exampleSummaryData.log2, method="quantile", transform="none")



###################################################
### code chunk number 17: normalise2
###################################################

exampleSummaryData.norm2 <- normaliseIllumina(channel(exampleSummaryData, "G.ul"), method="neqc", transform="none")



###################################################
### code chunk number 18: filter
###################################################

library(illuminaHumanv3.db)

ids <- as.character(featureNames(exampleSummaryData.norm))

qual <- unlist(mget(ids, illuminaHumanv3PROBEQUALITY, ifnotfound=NA))

table(qual)

rem <- qual == "No match" | qual == "Bad" | is.na(qual)

exampleSummaryData.filt <- exampleSummaryData.norm[!rem,]

dim(exampleSummaryData.filt)

boxplot(exampleSummaryData.norm, probeFactor = "PROBEQUALITY", sampleFactor="SampleFac") + opts(axis.text.x=theme_text(angle=45, hjust=1)) +  xlab("Probe Quality") + ylab("Log2 Intensity")



###################################################
### code chunk number 19: deanalysis
###################################################
library(limma)

rna <- factor(pData(exampleSummaryData)[,"SampleFac"])

design <- model.matrix(~0+rna)
colnames(design) <- levels(rna)
aw <- arrayWeights(exprs(exampleSummaryData.filt), design)
aw
fit <- lmFit(exprs(exampleSummaryData.filt), design, weights=aw)
contrasts <- makeContrasts(UHRR-Brain, levels=design)
contr.fit <- eBayes(contrasts.fit(fit, contrasts))
topTable(contr.fit, coef=1)



###################################################
### code chunk number 20: annotation
###################################################

ids2 <- featureNames(exampleSummaryData.filt)

chr <- mget(ids2, illuminaHumanv3CHR, ifnotfound = NA)
cytoband<- mget(ids2, illuminaHumanv3MAP, ifnotfound = NA)
refseq <- mget(ids2, illuminaHumanv3REFSEQ, ifnotfound = NA)
entrezid <- mget(ids2, illuminaHumanv3ENTREZID, ifnotfound = NA)
symbol <- mget(ids2, illuminaHumanv3SYMBOL, ifnotfound = NA)
genename <- mget(ids2, illuminaHumanv3GENENAME, ifnotfound = NA)

anno <- data.frame(Ill_ID = ids2, Chr = as.character(chr),
           Cytoband = as.character(cytoband), RefSeq = as.character(refseq),
           EntrezID = as.numeric(entrezid), Symbol = as.character(symbol),
           Name = as.character(genename))

contr.fit$genes <- anno
topTable(contr.fit)



###################################################
### code chunk number 21: readBeadSummary (eval = FALSE)
###################################################
## 
## library(beadarray)
## dataFile = "AsuragenMAQC-probe-raw.txt"
## qcFile = "AsuragenMAQC-controls.txt"
## BSData = readBeadSummaryData(dataFile = dataFile,
## qcFile = qcFile, controlID = "ProbeID",
## skip = 0, qc.skip = 0, qc.columns = list(exprs = "AVG_Signal",
## Detection = "Detection Pval"))
## 


###################################################
### code chunk number 22: options
###################################################
options(width = 80)


###################################################
### code chunk number 23: sessionInfo
###################################################
sessionInfo()