###################################################
### chunk number 1:  eval=FALSE
###################################################
## options(width=70)


###################################################
### chunk number 2: Import
###################################################
library(beadarray)
targets = read.table("targets.txt", sep="\t", header=TRUE, as.is=TRUE)
targets
beadInfo = read.csv("genes.csv", as.is=TRUE)
BLData = readIllumina(arrayNames=targets$arrayID, textType=".csv", targets=targets, backgroundMethod="none")


###################################################
### chunk number 3: BLData
###################################################
slotNames(BLData)

an=arrayNames(BLData)
an

numBeads(BLData)

pData(BLData)

colnames(BLData[[1]])

BLData[[1]]$G[1:5]
BLData[[1]]$Gb[1:5]

getArrayData(BLData, array = 1, which = "G", log = FALSE)[1:5]
getArrayData(BLData, array = 1, which = "Gb", log = FALSE)[1:5]


###################################################
### chunk number 4: Boxplots
###################################################
BLData.bc = backgroundCorrect(BLData)

negs = NULL
for(i in 1:10) {
  negs[i] = sum(getArrayData(BLData.bc, array=i, log=FALSE)<0)
}
negs

ylim = c(4,16)
par(mfrow=c(1,3), mai=c(1.5,0.6,0.2,0.1))
boxplotBeads(BLData,las=2, outline=FALSE, ylim=ylim, main="fg", ylab=expression(log[2](intensity)))
boxplotBeads(BLData,las=2, whatToPlot ="Gb", outline=FALSE, ylim=ylim, main="bg")
boxplotBeads(BLData.bc,las=2, whatToPlot ="G", outline=FALSE, ylim=ylim,  main="fg-bg")


###################################################
### chunk number 5: imageplots
###################################################
par(mfrow=c(2,5))
zlim = c(6,16)
for(i in 1:10){
  imageplot(BLData.bc, array=i, nrow=50, ncol=50, high="red", low="yellow", zlim=zlim, main=an[i])
}

x11()
par(mfrow=c(2,5))
zlim = c(-1,1)
for(i in 1:10){
  imageplot(BLData.bc, whatToPlot="residG", array=i, nrow=50, ncol=50, high="red", low="yellow", zlim=zlim, main=an[i])
}


###################################################
### chunk number 6: outliers
###################################################
par(mfrow=c(2,5))
for(i in 1:10){
  o=findAllOutliers(BLData.bc, array=i)
  plotBeadLocations(BLData.bc, array=i, BeadIDs=o, main=an[i], SAM=TRUE, pch=".")
}

outliers = NULL
for(i in 1:10) {
  outliers[i] = length(findAllOutliers(BLData.bc, array=i))
}
x11()
par(mai=c(2,1,0.2,0.1))
barplot(outliers/numBeads(BLData.bc)*100, main="Outliers per array", ylab="%", las=2, names=an)


###################################################
### chunk number 7: createBeadSummaryData
###################################################
BSData = createBeadSummaryData(BLData.bc, imagesPerArray=1)

slotNames(BSData)
names(assayData(BSData))

dim(exprs(BSData))
dim(se.exprs(BSData))

exprs(BSData)[1:10,1:2]
se.exprs(BSData)[1:10, 1:2]
pData(BSData)

par(mai=c(2,1,0.2,0.1), mfrow=c(1,2))
boxplot(as.data.frame(log2(exprs((BSData)))), ylab=expression(log[2](intensity)), las=2)
boxplot(as.data.frame(NoBeads(BSData)),  ylab="number of beads", ylim=c(0,60), las=2)


###################################################
### chunk number 8: matchIDs
###################################################
rownames(exprs(BSData))[1:5]
index = match(rownames(exprs(BSData)), beadInfo$ProbeID)
illuminaIDs = beadInfo$Target[index]
matching = !is.na(illuminaIDs)
rownames(exprs(BSData))[matching] = illuminaIDs[matching]
rownames(exprs(BSData))[1:5]


###################################################
### chunk number 9: plotMAXY
###################################################
plotMAXY(exprs(BSData), arrays=1:5)

x11()
plotMAXY(exprs(BSData), arrays=6:10)


###################################################
### chunk number 10: normalising
###################################################
BSData.log2.quantile = normaliseIllumina(BSData, method="quantile", transform="log2")

BSData.vst.quantile = normaliseIllumina(BSData, method="quantile", transform="vst")

plotMAXY(exprs(BSData.log2.quantile), arrays=1:5, log=FALSE)
plotMAXY(exprs(BSData.vst.quantile), arrays=1:5, log=FALSE)


###################################################
### chunk number 11: snpexample
###################################################
data(BLData)

numBeads(BLData)

colnames(BLData[[1]])

getArrayData(BLData, array = 1, which = "G", log = FALSE)[1:5]
getArrayData(BLData, array = 1, which = "R", log = FALSE)[1:5]


###################################################
### chunk number 12: snpqc
###################################################
x11()
par(mfrow=c(1,2))
plotRG(BLData, arrays=1)
plotRG(BLData, arrays=2)

x11()
plotBeadDensities(BLData, what=c("G", "R"), col=rep(c(4,7), times=4), lwd=2)
legend(13,1.1,legend=c("Cy3", "Cy5"), col=c(4,7), lwd=2)


###################################################
### chunk number 13: snpplots
###################################################
x11()
par(mfrow=c(1,2))
plotRG(BLData, ProbeIDs="1037", arrays=1:4, main="SNP 1037", pch=16)
plotRG(BLData, ProbeIDs="913", arrays=1:4, main="SNP 913", pch=16)


###################################################
### chunk number 14: snpsummary
###################################################
RG = createBeadSummaryData(BLData, what="RG", log=TRUE)
class(RG)
slotNames(RG)
names(RG@assayData)

A = createBeadSummaryData(BLData, what="A", log=TRUE)
M = createBeadSummaryData(BLData, what="M", log=TRUE)
dim(exprs(M))

x11()
boxplot(as.data.frame(exprs(M)), ylab=expression(log[2](R/G)))


###################################################
### chunk number 15: sessionInfo
###################################################
sessionInfo()