Entering edit mode
pingzhao Hu
▴
210
@pingzhao-hu-685
Last seen 10.2 years ago
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