##
## RMA-MS sample code using a dummy example
##  - this works as of BioC 2.3 in R-2.8.1 (Windows)
##

## First load data -- create an AffyBatch object called Data

library(ALLMLL); data(MLL.B)
Data <- MLL.B #AffyBatch object, 20 arrays



## Now define some functions (shouldn't need to change any of this)

######################################################################################
# RMA-MS: abatch is AffyBatch object, gn.sig is a vector of signal-possible geneNames
rma.ms <- function(abatch,gn.sig)
{ rma.abatch <- bg.correct(abatch, method="ms", gn.sig=gn.sig)
rma.abatch.norm <- qnorm.subset(rma.abatch,gn.sig)
out.eset <- expresso(rma.abatch.norm,bg.correct=FALSE,normalize=FALSE,
pmcorrect.method="pmonly",summary.method="medianpolish",summary.subset=gn.sig)
return(out.eset)
}
# RMA-MS background correction; based on bg.correct.rma
bg.correct.ms <- function (object, gn.sig, ...)
{ pn <- probeNames(object); t.sig <- is.element(pn,gn.sig)
pm.mat <- pm(object); pm.sig <- apply(pm.mat, 2, bg.adjust.ms, t.sig=t.sig)
pm.mat <- pm.sig; pm(object) <- pm.mat
return(object)
}
bg.adjust.ms <- function (pm, n.pts = 2^14, t.sig, ...)
{ param <- bg.parameters.ms(pm, n.pts, t.sig)
b <- param$sigma; pm <- pm - param$mu - param$alpha * b^2
pm + b * ((1/sqrt(2 * pi)) * exp((-1/2) * ((pm/b)^2)))/pnorm(pm/b)
}
bg.parameters.ms <- function (pm, n.pts = 2^14, t.sig)
{ pm.bg <- pm[!t.sig]; pm.sig <- pm[t.sig]
mu.bg <- median(pm.bg); sd.bg <- mad(pm.bg)
alpha.sig <- mean(pm.sig[pm.sig>mu.bg])
list(alpha = 1/alpha.sig, mu = mu.bg, sigma = sd.bg) # Note 1/alpha.sig
}
upDate.bgcorrect.methods(c(bgcorrect.methods(),"ms"))
# Quantile normalize a subset of genes; based on normalize.AffyBatch.quantiles
qnorm.subset <- function(abatch,subset)
{ pms <- unlist(pmindex(abatch,subset))
noNA <- rowSums(is.na(intensity(abatch)[pms, , drop = FALSE])) == 0
pms <- pms[noNA]
intensity(abatch)[pms, ] <- normalize.quantiles(intensity(abatch)[pms, , drop = FALSE],
copy = FALSE)
return(abatch)
}
######################################################################################



## Now define signal-possible geneNames (gn.sp) and call rma.ms function

library(affy); library(preprocessCore)
gn.sp <- geneNames(Data)[1:1000] # vector of signal-possible geneNames (probeset ids)
eset.ms <- rma.ms(Data,gn.sp) 

## This final eset.ms object is an ExpressionSet object, with only 
## expression estimates for signal-possible geneNames (in gn.sp)