This message is in MIME format. The first part should be readable
text,
while the remaining parts are likely unreadable without MIME-aware
tools.
Send mail to mime@docserver.cac.washington.edu for more info.
--Boundary_(ID_a6VjjKAmfoOJ3ucipeawFg)
Content-type: TEXT/PLAIN; charset=US-ASCII
Content-transfer-encoding: 7BIT
ive done 184. it used a lot of RAM though. how much RAM do you have?
for
the next version there will be a method that only uses PMs. this
reduces
the data by 50%. im attaching a function that might help you with a
work
around... if you can get more RAM you can try
1-reading in using read.celfile
2-keeping only the pms using 'pmormm'
3-using a version of the attached function to get expression.
hope this helps,
rafael
On Mon, 6 May 2002, Isaac Neuhaus wrote:
> Have anybody processed more than 64 cel files. I have an experiment
with
> 66 cel files and when it start processing file number 65 the
function
> ReadAffy dies. Apparently is dying in this line:
>
> aux <-
> as.vector(read.celfile(CELfiles[i],compress=compress.cel,sd=FALSE)$i
ntensity)
>
> I have tried some combinations of the 64 files just to make sure
that
> there is no corruption in one of the cel files and it works just
fine...
> up to to 64.
>
> Any suggestions?
>
> Isaac
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor@stat.math.ethz.ch
>
http://www.stat.math.ethz.ch/mailman/listinfo/bioconductor
>
--Boundary_(ID_a6VjjKAmfoOJ3ucipeawFg)
Content-id: <pine.lnx.4.33.0205061028390.27828@localhost.localdomain>
Content-type: TEXT/PLAIN; charset=US-ASCII; name=rma.R
Content-transfer-encoding: 7BIT
Content-disposition: attachment; filename=rma.R
Content-description:
##this function does the expression for our method: RMA. not flexible
rma <- function(object,subset=NULL, phenodata=NULL, annotation=NULL,
description=NULL, notes=NULL, verbose=TRUE, ...){
bg.parameters2 <- function(pm,n.pts=2^14){
max.density <- function(x,n.pts){
aux <- density(x,kernel="epanechnikov",n=n.pts)
aux$x[order(-aux$y)[1]]
}
pmbg <- max.density(pm,n.pts) ##Log helps detect mode
bg.data <- pm[pm < pmbg]
##do it again to really get the mode
pmbg <- max.density(bg.data,n.pts)
bg.data <- pm[pm < pmbg]
bg.data <- bg.data - pmbg
bgsd <- sqrt(sum(bg.data^2)/(length(bg.data)-1))*sqrt(2)#/.85
sig.data <- pm[pm > pmbg]
sig.data <- sig.data-pmbg
expmean <- max.density(sig.data,n.pts)
alpha <- 1/expmean
mubg <- pmbg
list(alpha=alpha,mu=mubg,sigma=bgsd)
}
bg.adjust2 <- function(pm,n.pts=2^14){
param <- bg.parameters2(pm,n.pts)
b <- param$sigma
a <- pm - param$mu - param$alpha*b^2
a + b*((1./sqrt(2*pi))*exp((-1./2.)*((a/b)^2)))/pnorm(a/b)
}
if(verbose) cat("Background correcting\n")
x <- apply(pm(object),2,bg.adjust2)
if(verbose) cat("Normalizing Data\n")
x <- normalize.quantiles(x)
if(verbose) cat("Preparing Data\n")
if(is.null(subset)){
data.lst <-split(as.vector(x),factor(rep(probeNames(object),nprobe
s(object)/dim(x)[1]*dim(x)[2])))
}
else{
if(is.numeric(subset)) subset <- geneNames(object)[subset]
Index <- which(probeNames(object)%in%subset)
x <- x[Index,]
Names <- probeNames(object)[Index]
data.lst <- split(as.vector(x),
factor(rep(Names,length(Names)/dim(x)[1]*dim(x)[2] )))
}
data.lst<-lapply(data.lst,matrix,ncol=nchips(object))
if(verbose) cat("Computing expression. This may take a while.\n")
e <- t(sapply(data.lst,medianpolish,...))
if(is.null(phenodata)) phenodata <- phenoData(object)
if(is.null(annotation)) annotation <- annotation(object)
if(is.null(description)) description <- description(object)
if(is.null(notes)) notes <- notes(object)
exprs <- e[,1:nchips(object)]
colnames(exprs) <- sampleNames(object)
se.exprs <- e[,(nchips(object)+1):(2*nchips(object))]
colnames(se.exprs) <- sampleNames(object)
new("exprSet",
exprs=exprs,
se.exprs=se.exprs,
phenoData=phenodata,
annotation=annotation,
description=description,
notes=notes)
}
--Boundary_(ID_a6VjjKAmfoOJ3ucipeawFg)--
