### R code from vignette source 'blima.rnw'
### Encoding: UTF-8

###################################################
### code chunk number 1: style-Sweave
###################################################
BiocStyle::latex(use.unsrturl=FALSE)


###################################################
### code chunk number 2: blima.rnw:21-26
###################################################
library(blima)
library(blimaTestingData)
data(blimatesting)
library(Biobase)
library(xtable)


###################################################
### code chunk number 3: blima.rnw:29-47
###################################################
array1stats = chipArrayStatistics(blimatesting[[1]], includeBeadStatistic=TRUE, 
	excludedOnSDMultiple=3)
array1pheno = pData(blimatesting[[1]]@experimentData$phenoData)
array1stats = data.frame(array1pheno$Name, array1stats)
colnames(array1stats)[1] <- "Array";
table = xtable(array1stats, align="|c|c|c|c|c|c|c|c|c|c|", caption="Array 1 statistic.")
digits(table)[c(2,3)]<-0
digits(table)[c(4:9)]<-1
print(table, include.rownames=FALSE)
array2stats = chipArrayStatistics(blimatesting[[2]], includeBeadStatistic=TRUE, 
	excludedOnSDMultiple=3)
array2pheno = pData(blimatesting[[2]]@experimentData$phenoData)
array2stats = data.frame(array2pheno$Name, array2stats)
colnames(array2stats)[1] <- "Array";
table = xtable(array2stats, align="|c|c|c|c|c|c|c|c|c|c|", caption="Array 2 statistic.")
digits(table)[c(2,3)]<-0
digits(table)[c(4:9)]<-1
print(table, include.rownames=FALSE)


###################################################
### code chunk number 4: blima.rnw:51-61
###################################################
library(illuminaHumanv4.db)
adrToIllumina = toTable(illuminaHumanv4ARRAYADDRESS)
adrToIllumina = adrToIllumina[, c("ArrayAddress", "IlluminaID")]
colnames(adrToIllumina) = c("Array_Address_Id", "Probe_Id")
illuminaToSymbol = toTable(illuminaHumanv4SYMBOLREANNOTATED)
adrToSymbol = merge(adrToIllumina, illuminaToSymbol, by.x="Probe_Id", by.y="IlluminaID")
adrToSymbol = adrToSymbol[,c("Array_Address_Id", "SymbolReannotated")]
colnames(adrToSymbol) = c("Array_Address_Id", "Symbol")
negIl = mappedLkeys(revmap(illuminaHumanv4REPORTERGROUPNAME)["negative"])
negAdr = mappedRkeys(illuminaHumanv4ARRAYADDRESS[negIl])


###################################################
### code chunk number 5: blima.rnw:64-72
###################################################
if(exists("annotationHumanHT12V4"))
{
    adrToIllumina = annotationHumanHT12V4$Probes[, c("Array_Address_Id", "Probe_Id")]
    adrToSymbol = annotationHumanHT12V4$Probes[, c("Array_Address_Id", "Symbol")]
    negAdr = unique(annotationHumanHT12V4$Controls[
       annotationHumanHT12V4$Controls$Reporter_Group_Name=="negative", 
       "Array_Address_Id"])
}


###################################################
### code chunk number 6: blima.rnw:75-79 (eval = FALSE)
###################################################
## blimatestingall = bacgroundCorrect(blimatesting, channelBackground = "GrnB", 
##         channelBackgroundFilter="bgf")
## blimatestingall = nonPositiveCorrect(blimatestingall, normalizationMod=NULL, 
##         channelCorrect="GrnF",  channelBackgroundFilter="bgf", channelAndVector="bgf")


###################################################
### code chunk number 7: blima.rnw:82-84 (eval = FALSE)
###################################################
## blimatestingall = backgroundChannelSubtract(blimatestingall, normalizationMod = NULL, 
##         channelSubtractFrom = "GrnF", channelSubtractWhat = "GrnB", channelResult = "BGS")


###################################################
### code chunk number 8: blima.rnw:87-90 (eval = FALSE)
###################################################
## blimatestingall = xieBacgroundCorrect(blimatestingall, normalizationMod = NULL, 
##         negativeArrayAddresses=negAdr, channelCorrect="GrnF", channelResult="XIE",
##         channelInclude="bgf")


###################################################
### code chunk number 9: blima.rnw:93-95 (eval = FALSE)
###################################################
## blimatestingall = varianceBeadStabilise(blimatestingall, quality="GrnF", 
##         channelInclude="bgf", channelOutput="vst")


###################################################
### code chunk number 10: blima.rnw:98-101 (eval = FALSE)
###################################################
## blimatestingall = selectedChannelTransform(blimatestingall, normalizationMod=NULL, 
##         channelTransformFrom="GrnF", channelResult="LOG", 
##         transformation=log2TransformPositive)


###################################################
### code chunk number 11: blima.rnw:104-106 (eval = FALSE)
###################################################
## blimatestingall = quantileNormalize(blimatestingall, normalizationMod=NULL, 
##         channelNormalize="vst", channelOutput="qua", channelInclude="bgf")


###################################################
### code chunk number 12: blima.rnw:110-125
###################################################
data("blimatesting")
groups1 = "A";
groups2 = "E";
sampleNames = list()
groups1Mod = list()
groups2Mod = list()
processingMod = list()
for(i in 1:length(blimatesting))
{
	p = pData(blimatesting[[i]]@experimentData$phenoData)
	groups1Mod[[i]] = p$Group %in% groups1;
	groups2Mod[[i]] = p$Group %in% groups2;
	processingMod[[i]] = p$Group %in% c(groups1, groups2);
	sampleNames[[i]] = p$Name
}


###################################################
### code chunk number 13: blima.rnw:128-136
###################################################
blimatesting = bacgroundCorrect(blimatesting, normalizationMod = processingMod, 
		channelBackgroundFilter="bgf")
blimatesting = nonPositiveCorrect(blimatesting, normalizationMod = processingMod, 
		channelCorrect="GrnF",  channelBackgroundFilter="bgf", channelAndVector="bgf")
blimatesting = varianceBeadStabilise(blimatesting, normalizationMod = processingMod, 
		quality="GrnF", channelInclude="bgf", channelOutput="vst")
blimatesting = quantileNormalize(blimatesting, normalizationMod = processingMod, 
		channelNormalize="vst", channelOutput="qua", channelInclude="bgf")


###################################################
### code chunk number 14: blima.rnw:139-168
###################################################
probeTest <- doProbeTTests(blimatesting, groups1Mod, groups2Mod, 
		transformation=NULL, quality="qua", channelInclude="bgf")
beadTest <- doTTests(blimatesting, groups1Mod, groups2Mod, 
		transformation=NULL, quality="qua", channelInclude="bgf")
probeTestID = probeTest[,"ProbeID"]
beadTestID = beadTest[,"ProbeID"]
probeTestFC = abs(probeTest[,"mean1"]-probeTest[,"mean2"])
beadTestFC = abs(beadTest[,"mean1"]-beadTest[,"mean2"])
probeTestP = probeTest[,"adjustedp"]
beadTestP = beadTest[,"adjustedp"]
probeTestMeasure = (1-probeTestP)*probeTestFC
beadTestMeasure = (1-beadTestP)*beadTestFC
probeTest = cbind(probeTestID, probeTestMeasure)
beadTest = cbind(beadTestID, beadTestMeasure)
colnames(probeTest) <- c("ArrayAddressID", "difexPL")
colnames(beadTest) <- c("ArrayAddressID", "difexBL")
tocmp <- merge(probeTest, beadTest)
tocmp = merge(tocmp, adrToSymbol, by.x="ArrayAddressID", by.y="Array_Address_Id")
tocmp = tocmp[, c("ArrayAddressID", "Symbol", "difexPL", "difexBL")]
sortPL = sort(-tocmp[,"difexPL"], index.return=TRUE)$ix
sortBL = sort(-tocmp[,"difexBL"], index.return=TRUE)$ix
beadTop10 = tocmp[sortBL[1:10],]
probeTop10 = tocmp[sortPL[1:10],]
beadTop10 = xtable(beadTop10, align="|c|c|c|c|c|", caption="Top 10 probes on bead level.")
probeTop10 = xtable(probeTop10, align="|c|c|c|c|c|", caption="Top 10 probes on probe level.")
digits(beadTop10)[2] = 0
digits(probeTop10)[2] = 0
print(beadTop10, include.rownames=FALSE)
print(probeTop10, include.rownames=FALSE)


###################################################
### code chunk number 15: blima.rnw:174-201
###################################################
nonnormalized = createSummarizedMatrix(blimatesting, spotsToProcess=processingMod, quality="GrnF", channelInclude="bgf", 
        annotationTag="Name")
nonnormalized = merge(nonnormalized, adrToIllumina, by.x="ProbeID", by.y="Array_Address_Id")
nonnormalized = nonnormalized[, c(10, 2:9)]
colnames(nonnormalized)[1] = "ID_REF"
for(i in 2:9)
{
    colnames(nonnormalized)[i] = sprintf("%s", colnames(nonnormalized)[i]) 
}
table = head(nonnormalized)
table = xtable(table, align="|c|c|c|c|c|c|c|c|c|c|", caption="Head of nonnormalized data.")
digits(table)[c(2:9)]<-1
print(table, include.rownames=FALSE)

normalized = createSummarizedMatrix(blimatesting, spotsToProcess=processingMod, quality="qua", channelInclude="bgf", 
        annotationTag="Name")
normalized = merge(normalized, adrToIllumina, by.x="ProbeID", by.y="Array_Address_Id")
normalized = normalized[, c(10, 2:9)]
colnames(normalized)[1] = "ID_REF"
for(i in 2:9)
{
    colnames(normalized)[i] = sprintf("%s", colnames(normalized)[i])
}
table = head(normalized)
table = xtable(table, align="|c|c|c|c|c|c|c|c|c|c|", caption="Head of normalized data.")
digits(table)[c(2:10)]<-3
print(table, include.rownames=FALSE)