[BioC] Max Cel Files

[Rafael A. Irizarry]

  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)$intensity)
> 
> 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),nprobes(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)--