affyContam.R
In affyContam: structured corruption of affymetrix cel file data

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

###################################################
### code chunk number 1: lkd
###################################################
library(affydata)
data(Dilution)
image(Dilution[,1])


###################################################
### code chunk number 2: dodef
###################################################
library(affyContam)
dilc = setCircRegion(Dilution, chip=1)


###################################################
### code chunk number 3: lkp (eval = FALSE)
###################################################
## image(dilc[,1])


###################################################
### code chunk number 4: getl (eval = FALSE)
###################################################
## library(affyMvout)
## library(affy)
## library(SpikeIn)
## data(SpikeIn133)
## library(mdqc)
## library(affyContam)
## library(limma)


###################################################
### code chunk number 5: doperm (eval = FALSE)
###################################################
## s12 = SpikeIn133[,1:12]
## s12rma = rma(s12)
## mads = apply(exprs(s12rma),1,mad)
## kp = which(mads > quantile(mads,.95))
## kppn = featureNames(s12rma)[kp]
## # these are the 18 genes found to be mostly monotone over 12 chips
## mostmr = c("203508_at", "204563_at", "204513_s_at", "204205_at", "204959_at", 
## "207655_s_at", "204836_at", "205291_at", "209795_at", "207777_s_at", 
## "204912_at", "205569_at", "207160_at", "205692_s_at", "212827_at", 
## "AFFX-LysX-3_at", "AFFX-PheX-3_at", "AFFX-ThrX-3_at")
## 


###################################################
### code chunk number 6: cod (eval = FALSE)
###################################################
## fullrun = function( abatch, arma, contFun, filtpn, targpn, chips=1, ... ) {
## # assess detectability in original data
##  dvec = (1:ncol(exprs(abatch)))
##  des = model.matrix(~dvec)
##  af1 = lmFit( arma[filtpn,], des, method="robust", maxit=300 )
##  eaf1 = eBayes(af1)
##  orig.tt = eaf1$t[targpn,2]
## # contaminate
##  cbat = contFun(abatch, chip=chips[1], ...)
##  if (length(chips)>1) {
##    for (i in 2:(length(chips)))
##      cbat = contFun(cbat, chip=chips[i], ...)
##    }
## # assess detectability in contaminated data
##  crma = rma(cbat)[filtpn,]
##  dvec = (1:ncol(exprs(abatch)))
##  des = model.matrix(~dvec)
##  cf1 = lmFit( crma, des, method="robust", maxit=300 )
##  ecf1 = eBayes(cf1)
##  contam.tt = ecf1$t[targpn,2]
## # now test for outliers
##  caos = ArrayOutliers(cbat)
##  if (nrow(caos[[1]]) < 1) {
##     warning("no outliers by affyMvout")
##     return(list(aos=caos, md=mdqc(caos[[3]][,2:10])))
##     }
##  todrop = as.numeric(rownames(caos[[1]]))
##  cbatf = cbat[,-todrop]
## # assess detectability in repaired data
##  frma = rma(cbatf)[filtpn,]
##  dvec = (1:ncol(exprs(abatch)))[-todrop]
##  des = model.matrix(~dvec)
##  f1 = lmFit( frma, des, method="robust", maxit=300 )
##  ef1 = eBayes(f1)
##  repair.tt = ef1$t[targpn,2]
## # compute the mdqc result
##  md = mdqc(caos[[3]][,2:10], robust="MCD")
##  list(orig=orig.tt, contam=contam.tt, repair=repair.tt, md=md, todrop=todrop)
## }
##  


###################################################
### code chunk number 7: ddd (eval = FALSE)
###################################################
## ff = fullrun( s12, s12rma, setCircRegion, kppn, mostmr, chips=1:2 )


###################################################
### code chunk number 8: scr (eval = FALSE)
###################################################
## scr.40 = function (x, chip = 1, center = c(150, 150), rad = 75, vals = 30, 
##     valgen = NULL) 
## {
##     cdfname = paste(annotation(x), "cdf", sep = "")
##     require(cdfname, character.only = TRUE, quietly = TRUE)
##     xext = seq(center[1] - rad, center[1] + rad)
##     yext = seq(center[2] - rad, center[2] + rad)
##     badco = expand.grid(xext, yext)
##     badco = badco[(badco[, 1] - center[1])^2 + (badco[, 2] - 
##         center[2])^2 < rad^2, ]
##     indsbad = apply(badco, 1, function(x) xy2indices(x[1], x[2], 
##         cdf = cdfname))
##     if (is.null(valgen)) 
##         exprs(x)[indsbad, chip] = vals
##     else exprs(x)[indsbad, chip] = valgen(length(indsbad))
##     x
## }


###################################################
### code chunk number 9: scr2 (eval = FALSE)
###################################################
## scr.20k = function (x, chip = 1, center = c(500, 500), rad = 75, vals = 30000,
##     valgen = NULL)
## {
##     cdfname = paste(annotation(x), "cdf", sep = "")
##     require(cdfname, character.only = TRUE, quietly = TRUE)
##     xext = seq(center[1] - rad, center[1] + rad)
##     yext = seq(center[2] - rad, center[2] + rad)
##     badco = expand.grid(xext, yext)
##     badco = badco[(badco[, 1] - center[1])^2 + (badco[, 2] -
##         center[2])^2 < rad^2, ]
##     indsbad = apply(badco, 1, function(x) xy2indices(x[1], x[2],
##         cdf = cdfname))
##     if (is.null(valgen))
##         exprs(x)[indsbad, chip] = vals
##     else exprs(x)[indsbad, chip] = valgen(length(indsbad))
##     x
## }


###################################################
### code chunk number 10: dov (eval = FALSE)
###################################################
## incvarCircRegion = function(x, chip=1, center=c(150,500), rad=100, fac=3) {
##  tmp = fac*getCircRegion(x, chip, center, rad)
##  setCircRegion(x, chip, center, rad, vals=tmp)
## }


###################################################
### code chunk number 11: dov (eval = FALSE)
###################################################
## incvarRectRegion = function(x, chip=1, xinds=350:700, yinds=1:700, fac=3) {
##  tmp = fac*getRectRegion(x, chip, xinds, yinds)
##  setRectRegion(x, chip, xinds, yinds, vals=tmp)
## }


###################################################
### code chunk number 12: compo (eval = FALSE)
###################################################
## tryout = scr.40(s12)
## tryout = scr.20k(tryout)
## tryout = incvarCircRegion(tryout)
## fin = incvarRectRegion(tryout)
## png(file="lkcomp.png")
## image(fin[,1], main="composite contamination")
## dev.off()


###################################################
### code chunk number 13: doruns (eval = FALSE)
###################################################
## d1b = fullrun( s12, s12rma, scr.40, kppn, mostmr, chips=1 )
## save(d1b, file="d1b.rda")
## d2b = fullrun( s12, s12rma, scr.40, kppn, mostmr, chips=1:2 )
## save(d2b, file="d2b.rda")
## d3b = fullrun( s12, s12rma, scr.40, kppn, mostmr, chips=c(1:2,11) )
## save(d3b, file="d3b.rda")
## H1b = fullrun( s12, s12rma, scr.20k, kppn, mostmr, chips=1 )
## save(H1b, file="H1b.rda")
## H2b = fullrun( s12, s12rma, scr.20k, kppn, mostmr, chips=1:2 )
## save(H2b, file="H2b.rda")
## H3b = fullrun( s12, s12rma, scr.20k, kppn, mostmr, chips=c(1:2,11) )
## save(H3b, file="H3b.rda")
## I1b = fullrun( s12, s12rma, incvarCircRegion, kppn, mostmr, chips=1 )
## save(I1b, file="I1b.rda")
## I2b = fullrun( s12, s12rma, incvarCircRegion, kppn, mostmr, chips=1:2 )
## save(I2b, file="I2b.rda")
## I3b = fullrun( s12, s12rma, incvarCircRegion, kppn, mostmr, chips=c(1:2,11) )
## save(I3b, file="I3b.rda")
## R1b = fullrun( s12, s12rma, incvarRectRegion, kppn, mostmr, chips=1 )
## save(R1b, file="R1b.rda")
## R2b = fullrun( s12, s12rma, incvarRectRegion, kppn, mostmr, chips=1:2 )
## save(R2b, file="R2b.rda")
## R3b = fullrun( s12, s12rma, incvarRectRegion, kppn, mostmr, chips=c(1:2,11) )
## save(R3b, file="R3b.rda")