inst/doc/ARRmNormalization.R

In ARRmNormalization: Adaptive Robust Regression normalization for Illumina methylation data

### R code from vignette source 'ARRmNormalization.Rnw'

###################################################
### code chunk number 1: dataLoading
###################################################
library(ARRmData)
data(greenControlMatrix)
data(redControlMatrix)
data(betaMatrix)
data(sampleNames)


###################################################
### code chunk number 2: printBetaMatrix
###################################################
betaMatrix[1:5,1:3]


###################################################
### code chunk number 3: printControlMatrices
###################################################
greenControlMatrix[1:5,1:3]
redControlMatrix[1:5,1:3]


###################################################
### code chunk number 4: probesLoading
###################################################
library(ARRmNormalization)
data(ProbesType)
head(ProbesType)


###################################################
### code chunk number 5: printSampleNames
###################################################
head(sampleNames)


###################################################
### code chunk number 6: packageLoading
###################################################
library(ARRmNormalization)


###################################################
### code chunk number 7: dataPackageLoading
###################################################
library(ARRmData)


###################################################
### code chunk number 8: dataLoading2
###################################################
data(betaMatrix)
data(greenControlMatrix)
data(redControlMatrix)
data(sampleNames)


###################################################
### code chunk number 9: getBackground
###################################################
backgroundInfo <- getBackground(greenControlMatrix, redControlMatrix)


###################################################
### code chunk number 10: printBackground
###################################################
head(backgroundInfo)


###################################################
### code chunk number 11: printSampleNames
###################################################
head(sampleNames)


###################################################
### code chunk number 12: getDesignInfo
###################################################
designInfo <- getDesignInfo(sampleNames)
head(designInfo)


###################################################
### code chunk number 13: explicitDesign
###################################################
myChipVector <- c(rep(1,12), rep(2, 12), rep(3, 12))
myPositionVector <- rep(seq(1:12), 3)


###################################################
### code chunk number 14: getExplicitDesign
###################################################
designInfo <- getDesignInfo(sampleNames=NULL, 
                  positionVector=myPositionVector, chipVector=myChipVector)
head(designInfo)


###################################################
### code chunk number 15: firstNorm (eval = FALSE)
###################################################
## normMatrix <- 
##     normalizeARRm(betaMatrix=betaMatrix, designInfo=designInfo,
##         backgroundInfo=backgroundInfo, outliers.perc=0.1)


###################################################
### code chunk number 16: goodProbes
###################################################
data(ProbesType) 
goodProbes <- as.character(unlist(ProbesType$Probe_Name))[seq(1,485577,2)]


###################################################
### code chunk number 17: normGoodProbes (eval = FALSE)
###################################################
## normMatrix <- 
##     normalizeARRm(betaMatrix=betaMatrix, designInfo=designInfo,
##         backgroundInfo=backgroundInfo, outliers.perc=0, goodProbes=goodProbes)


###################################################
### code chunk number 18: normWOChip (eval = FALSE)
###################################################
## normMatrix <- 
##     normalizeARRm(betaMatrix=betaMatrix, designInfo=designInfo,
##         backgroundInfo=backgroundInfo, outliers.perc=0, chipCorrection=F)


###################################################
### code chunk number 19: getQuantiles
###################################################
quantiles=getQuantiles(betaMatrix)
attributes(quantiles)


###################################################
### code chunk number 20: printQuantiles
###################################################
quantiles=getQuantiles(betaMatrix, goodProbes=goodProbes)
attributes(quantiles)


###################################################
### code chunk number 21: printQuantiles2
###################################################
quantiles$II[1:5,1:4]


###################################################
### code chunk number 22: quantilesPlots
###################################################
quantilePlots(quantiles, backgroundInfo, designInfo)


###################################################
### code chunk number 23: plot1
###################################################
myColors=rep(c("black","red"),10)
	cex=1
	pch=rep(c(20,18),10)
	lwd=1.5
	percentilesI=seq(1,100,5)
	percentilesII=seq(1,100,10)
	y<-seq(from=-0.1,to=1,by=11/100)
	quantilesList=list(quantiles$green, quantiles$red, quantiles$II)
backgroundList = 
  list(log(backgroundInfo$green), log(backgroundInfo$red), log(backgroundInfo$red*backgroundInfo$green))
	titles=c("Background effects on different percentiles - Inf I Green",
           "Background effects on different percentiles - Inf I Red",
           "Background effects on different percentiles - Infinium II")
j=3
if (j==1 || j==2){percentiles=percentilesI} else {percentiles=percentilesII}
		currentQuantiles=quantilesList[[j]]
    	background=backgroundList[[j]]
		x<-seq(from=min(background), 
           to=max(background)+0.5, 
               by=(max(background+0.5)-min(background))/10)


###################################################
### code chunk number 24: fig1plot
###################################################
plot(x,y,type="n",xlab="Log Background intensity",ylab="Beta value",main=titles[j])
		for (i in 1:length(percentiles)){
			currentP=percentiles[i]
			points(unlist(background), 
                 unlist(currentQuantiles[currentP,]),
                     col=myColors[i],
                         pch=pch[i],
                             cex=cex)
			abline(lm(unlist(currentQuantiles[currentP,])~background),col="black",lwd=lwd)
		}	


###################################################
### code chunk number 25: fig1
###################################################
plot(x,y,type="n",xlab="Log Background intensity",ylab="Beta value",main=titles[j])
		for (i in 1:length(percentiles)){
			currentP=percentiles[i]
			points(unlist(background), 
                 unlist(currentQuantiles[currentP,]),
                     col=myColors[i],
                         pch=pch[i],
                             cex=cex)
			abline(lm(unlist(currentQuantiles[currentP,])~background),col="black",lwd=lwd)
		}	


###################################################
### code chunk number 26: plotDyeBias
###################################################
myColors=rep(c("black","red"),10)
	cex=1
	pch=rep(c(20,18),10)
	lwd=1.5
	percentilesI=seq(1,100,5)
	percentilesII=seq(1,100,10)
	y<-seq(from=-0.1,to=1,by=11/100)
	quantilesList=list(quantiles$green,quantiles$red,quantiles$II)

	## Dye bias plot
	dyebias=log(backgroundInfo$green/backgroundInfo$red)
	currentQuantile=quantilesList[[3]]
	percentiles=percentilesII
	x <- seq(from=min(dyebias),
           to=max(dyebias)+0.5,
               by=(max(dyebias+0.5)-min(dyebias))/10)


###################################################
### code chunk number 27: fig2plot
###################################################
plot(x,y,type="n",xlab="Dye bias",ylab="Beta value",main="Dye bias effects on different percentiles - Infinium II" )
		for (i in 1:length(percentiles)){
			currentP=percentiles[i]
			points(unlist(dyebias), 
             unlist(currentQuantile[currentP,]), 
                 col=myColors[i],pch=pch[i],cex=cex)
			abline(lm(unlist(currentQuantile[currentP,])~dyebias),col="black",lwd=lwd)
		}	


###################################################
### code chunk number 28: fig2
###################################################
plot(x,y,type="n",xlab="Dye bias",ylab="Beta value",main="Dye bias effects on different percentiles - Infinium II" )
		for (i in 1:length(percentiles)){
			currentP=percentiles[i]
			points(unlist(dyebias), 
             unlist(currentQuantile[currentP,]), 
                 col=myColors[i],pch=pch[i],cex=cex)
			abline(lm(unlist(currentQuantile[currentP,])~dyebias),col="black",lwd=lwd)
		}	


###################################################
### code chunk number 29: positionPlots
###################################################
positionPlots(quantiles, designInfo, percentiles=c(25, 50, 75))


###################################################
### code chunk number 30: examplePositionPlots
###################################################
	quantilesList= list(quantiles$green, quantiles$red, quantiles$II)
	percentiles=c(75)
	#### To look at the position deviations from the mean
	grandMeansList<-list()
	centeredQuantiles<-list()
	for (i in 1:3){grandMeansList[[i]]=apply(quantiles[[i]],1,mean)}
	for (i in 1:3){centeredQuantiles[[i]]=quantiles[[i]]-grandMeansList[[i]]}
	chipCenteredQuantiles=centeredQuantiles
	chipInfo=as.numeric(designInfo$chipInfo)
	chipIndices=unique(as.numeric(designInfo$chipInfo))
	for (i in 1:3){
		currentQuantiles=chipCenteredQuantiles[[i]]
		for (j in 1:length(chipInfo)){
			chipQuantiles=currentQuantiles[,which(chipInfo==chipIndices[j])]
			mean=apply(chipQuantiles,1,mean)
			chipQuantiles=chipQuantiles-mean
			currentQuantiles[,which(chipInfo==chipIndices[j])]=chipQuantiles
		}
		chipCenteredQuantiles[[i]]=currentQuantiles
	}

	### Position lines
	order=order(designInfo$positionInfo)
	for (i in 1:3){
    chipCenteredQuantiles[[i]]= chipCenteredQuantiles[[i]][,order]
	}
	newpositions= designInfo$positionInfo[order]
	lines=c()
	for (i in 1:(length(newpositions)-1)){
	if (newpositions[i+1]!=newpositions[i]){lines=c(lines,i+0.5)}
	}

	## Tick marks for position axis
	ticks=c()
	tick1=lines[1]/2
	lastTick=(nrow(designInfo)+lines[11])/2
	ticks=tick1
	for (i in 2:11){ticks=c(ticks,(lines[i]+lines[i-1])/2)}
	ticks=c(ticks,lastTick)


###################################################
### code chunk number 31: fig3plot
###################################################
	for (percIndex in 1:length(percentiles)){
		currentPerc=percentiles[percIndex]
		percStr=paste("Positions Variations - ",currentPerc,"th Percentile")
		titles=c(paste(percStr,"Inf I Green"), 
                 paste(percStr,"Inf I Red"), 
                     paste(percStr,"Inf II"))
    
		for (k in 3:3){
			plot(chipCenteredQuantiles[[k]][currentPerc,], 
           xaxt='n',ylab="Beta value deviations",xlab="Chip Positions",
               cex=1,pch=20,main=titles[k])
    
			axis (1, at = ticks,labels=1:12)
			abline(h=0)
			for (j in 1:length(lines)){abline(v=lines[j],lty=2)}
		}
	}


###################################################
### code chunk number 32: fig3
###################################################
	for (percIndex in 1:length(percentiles)){
		currentPerc=percentiles[percIndex]
		percStr=paste("Positions Variations - ",currentPerc,"th Percentile")
		titles=c(paste(percStr,"Inf I Green"), 
                 paste(percStr,"Inf I Red"), 
                     paste(percStr,"Inf II"))
    
		for (k in 3:3){
			plot(chipCenteredQuantiles[[k]][currentPerc,], 
           xaxt='n',ylab="Beta value deviations",xlab="Chip Positions",
               cex=1,pch=20,main=titles[k])
    
			axis (1, at = ticks,labels=1:12)
			abline(h=0)
			for (j in 1:length(lines)){abline(v=lines[j],lty=2)}
		}
	}


###################################################
### code chunk number 33: getCoefficients
###################################################
coefficients <- getCoefficients(quantiles, designInfo, backgroundInfo, outliers.perc=0.02)
attributes(coefficients)


###################################################
### code chunk number 34: printCoefficients
###################################################
attributes(coefficients$II)