###################################################
### chunk number 1: options
###################################################

options(width = 70)



###################################################
### chunk number 2: prelim
###################################################
library(beadarray)
data(BLData)


###################################################
### chunk number 3: readIllumina eval=FALSE
###################################################
## 
## BLData = readIllumina(useImages=FALSE, illuminaAnnotation = "Humanv3")


###################################################
### chunk number 4: BLData
###################################################

data(BLData)


is(BLData)
class(BLData)
slotNames(BLData)

##Get the beadData for array-section 1

BLData[[1]][1:10,]

##Alternative using accessor function

getBeadData(BLData, array=1, what="Grn")[1:10]

##Get unique ProbeIDs. These are the ArrayAddressIDs

uIDs = unique(getBeadData(BLData, array=1, what="ProbeID"))

uIDs[1:10]



###################################################
### chunk number 5: transform
###################################################

log2(BLData[[1]][1:10,2])

logGreenChannelTransform

logGreenChannelTransform(BLData, array=1)[1:10]

logRedChannelTransform



###################################################
### chunk number 6: backgroundCorrection
###################################################

for(i in 1:10){
 	BLData = backgroundCorrectSingleSection(BLData, array=i)
}

head(BLData[[1]])




###################################################
### chunk number 7: BeadLevelBoxplots
###################################################

getBackgroundCorrectionIntensities = function(BLData, array){

	log2(getBeadData(BLData, array=array, what="Grn.bc"))

}

par(mfrow=c(1,2))
boxplot(BLData,las=2, outline=FALSE, ylim=c(4,12), main="Green Foreground")
boxplot(BLData,las=2, transFun =getBackgroundCorrectionIntensities, outline=FALSE, ylim=c(4,12), main="Green Foreground")



###################################################
### chunk number 8: Imageplots
###################################################

par(mfrow=c(1,3), mar = c(2,2,2,2))

imageplot(BLData, array=1, low="lightgreen", high="darkgreen", horizontal = FALSE ,squareSize=25)

imageplot(BLData, array=2, low="lightgreen", high="darkgreen", horizontal = FALSE ,squareSize=25)

imageplot(BLData, array=3, low="lightgreen", high="darkgreen", horizontal = FALSE ,squareSize=25)




###################################################
### chunk number 9: BASH
###################################################


bsh = BASH(BLData, array=1:10)



###################################################
### chunk number 10: bashResults
###################################################

bsh$QC

par(mfrow=c(1,2))
barplot(bsh$QC[,1], main="Number of beads masked")

barplot(bsh$QC[,2], main="Extended Score")



###################################################
### chunk number 11: savingBASH
###################################################

for(i in 1:10){

	BLData = setWeights(BLData, wts=bsh$wts[[i]], array=i)

}

BLData = insertSectionData(BLData, what="BASHQC", data = bsh$QC)



###################################################
### chunk number 12: Controls
###################################################

data(controlProfile)

head(controlProfile)
table(controlProfile[,2])



###################################################
### chunk number 13: positiveControls
###################################################

par(mfrow=c(2,5), mar = c(2,2,2,2))

for(i in 1:10){
	poscontPlot(BLData, array=i, controlProfile=controlProfile, ylim = c(10,16))
}


###################################################
### chunk number 14: positiveNegativeControls
###################################################
par(mfrow=c(2,5), mar = c(2,2,2,2))

for(i in 1:10){
	poscontPlot(BLData, array=i, controlProfile=controlProfile,positiveControlTags = c("housekeeping", "negative"),  ylim = c(10,16))
}


###################################################
### chunk number 15: qaSummary
###################################################

quickSummary(BLData, array=1, reporterIDs=controlProfile[,1], reporterTags=controlProfile[,2])


qcReport = makeQCTable(BLData, controlProfile=controlProfile)

head(qcReport)[,1:5]

BLData = insertSectionData(BLData, what="BeadLevelQC", data=qcReport)

names(BLData@sectionData)

for(i in 1:10){
print(controlProbeDetection(BLData, array = i, controlProfile=controlProfile, tagsToDetect=c("housekeeping", "biotin"), negativeTag="negative"))
}




###################################################
### chunk number 16: qcPipeline eval=FALSE
###################################################
## 
## expressionQCPipeline(BLData, qcDir="QC")
## 
## expressionQCPipeline(BLData, controlProfile=controlProfile, positiveControlTags=c("housekeeping", "biotin"), hybridisationTags = c("cy3_hyb", "high_stringency_hyb"), zlim=c(10,12), horizontal=F, negativeTag = "negative", tagsToDetect=list(Housekeeping = "housekeeping", Biotin = "biotin", hybridisation = "cy3_hyb"))
## 
## 


###################################################
### chunk number 17: outlierPlots
###################################################

par(mfrow=c(1,3), mar = c(2,2,2,2))

outlierplot(BLData, array=1, horizontal = FALSE)

outlierplot(BLData, array=2,  horizontal = FALSE)

outlierplot(BLData, array=3, horizontal = FALSE)



###################################################
### chunk number 18: createBeadSummaryData
###################################################

myMean = function(x) mean(x,na.rm=TRUE)
mySd = function(x) sd(x,na.rm=TRUE)

greenChannel = new("illuminaChannel", logGreenChannelTransform, illuminaOutlierMethod, myMean, mySd,"G")

BSData <- summarize(BLData, list(greenChannel))


myMedian = function(x) median(x, na.rm=TRUE)
mySe = function(x) sd(x, na.rm=TRUE)/sqrt(length(x))


greenChannel2 = new("illuminaChannel", logGreenChannelTransform, illuminaOutlierMethod, myMedian, mySe,"G")

BSData2 <- summarize(BLData, list(greenChannel2))



###################################################
### chunk number 19: showBSData
###################################################
BSData

head(exprs(BSData)[,1:4])
head(se.exprs(BSData)[,1:4])

BSData2

head(exprs(BSData2)[,1:4])
head(se.exprs(BSData2)[,1:4])





###################################################
### chunk number 20: multipleChannel
###################################################

greenBackgroundCorrected = new("illuminaChannel", getBackgroundCorrectionIntensities, illuminaOutlierMethod, myMean, mySd, "G.bc")

BSData.multChannel = summarize(BLData, channelList = list(greenChannel, greenBackgroundCorrected))

channelNames(BSData.multChannel)

G = channel(BSData.multChannel, "G")

G.bc = channel(BSData.multChannel, "G.bc")




###################################################
### chunk number 21: calculateDetection
###################################################


status = rep("regular", as.numeric(dim(BSData.multChannel)[1]))

negIDs = controlProfile[which(controlProfile[,2] == "negative"),1]

status[match(negIDs, featureNames(BSData.multChannel))] = "negative"

det = calculateDetection(G, status=status)

head(det)

Detection(G) = det




###################################################
### chunk number 22: sessionInfo
###################################################

sessionInfo()