inst/doc/blima.R

In blima: Tools for the preprocessing and analysis of the Illumina microarrays on the detector (bead) level

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

## --------------------------------------------------------------------------
library(blima)
library(blimaTestingData)
data(blimatesting)
library(Biobase)
library(xtable)

## ----results='asis'--------------------------------------------------------
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)

## --------------------------------------------------------------------------
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])

## --------------------------------------------------------------------------
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"])
}

## ----eval=F----------------------------------------------------------------
#  blimatestingall = bacgroundCorrect(blimatesting, channelBackground = "GrnB",
#          channelBackgroundFilter="bgf")
#  blimatestingall = nonPositiveCorrect(blimatestingall, normalizationMod=NULL,
#          channelCorrect="GrnF",  channelBackgroundFilter="bgf", channelAndVector="bgf")

## ----eval=F----------------------------------------------------------------
#  blimatestingall = backgroundChannelSubtract(blimatestingall, normalizationMod = NULL,
#          channelSubtractFrom = "GrnF", channelSubtractWhat = "GrnB", channelResult = "BGS")

## ----eval=F----------------------------------------------------------------
#  blimatestingall = xieBacgroundCorrect(blimatestingall, normalizationMod = NULL,
#          negativeArrayAddresses=negAdr, channelCorrect="GrnF", channelResult="XIE",
#          channelInclude="bgf")

## ----eval=F----------------------------------------------------------------
#  blimatestingall = varianceBeadStabilise(blimatestingall, quality="GrnF",
#          channelInclude="bgf", channelOutput="vst")

## ----eval=F----------------------------------------------------------------
#  blimatestingall = selectedChannelTransform(blimatestingall, normalizationMod=NULL,
#          channelTransformFrom="GrnF", channelResult="LOG",
#          transformation=log2TransformPositive)

## ----eval=F----------------------------------------------------------------
#  blimatestingall = quantileNormalize(blimatestingall, normalizationMod=NULL,
#          channelNormalize="vst", channelOutput="qua", channelInclude="bgf")

## --------------------------------------------------------------------------
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
}

## --------------------------------------------------------------------------
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")

## ----results='asis'--------------------------------------------------------
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)

## ----results='asis'--------------------------------------------------------
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)