[BioC] RMA

[pingzhao Hu]

I simulated PM and MM matrix (named pm.raw and mm.raw)
(say 100 * 16 *15 --100 genes, 16 probes and 15 chips)
Now I try to estimate expression values using RMA.

Essentially, RMA has 3 steps:

1) background correction:
pm.bgc<- apply(pm.raw, 2, bg.adjust)
this is not working.

------------------
 >Error in density(x, kernel = "epanechnikov", n = n.pts) :
----------------------------------------
So I use:


mypm<-bg.adjust(pm.raw,mm.raw)
pm.bgc<-mypm$pm

bg.adjust<-
function (pm, mm,n.pts = 2^14, ...)
{
     param <- bg.parameters(pm, mm,n.pts)
     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)
     list(pm=pm)
}

### patch suggested by Laurent Gautier 10/17/2002 ####

bg.parameters <- function(pm, mm, n.pts=2^14){

   max.density <- function(x,n.pts){
     aux <- density(x, kernel="epanechnikov", n=n.pts, na.rm=TRUE)
     aux$x[order(-aux$y)[1]]
   }
     mmbg <- max.density(mm,n.pts)
     pmbg <- max.density(pm,n.pts)
		
     bg.data <- mm[mm < mmbg]
     bg.data <- bg.data-mmbg
     bgsd <- sqrt(sum(bg.data^2)/(length(bg.data)-1))*sqrt(2)/.85
			
     sig.data <- pm[pm > pmbg]
     sig.data <- sig.data-pmbg
			
     alpha <- 1/mean(sig.data)
     mubg <- mmbg
     list(alpha=alpha,mu=mubg,sigma=bgsd)
}


Is this correct for background correction????



(2) normalization

pm.norm<-normalize.quantiles(pm.bgc)


(3) probe specific correction and expression summary

It seems that express.summary.stat can do this.
However, the input into this function is a ProbeSet object.
Now I have pm.norm and mm.raw, how can I create this ProbeSet object????

or how can I get summarized expression values from normalized probel data???

Thanks

Pingzhao