[BioC] Farms/Call.INI using oligo package?

Djork-Arné Clevert okko at clevert.de
Fri Nov 15 15:46:10 CET 2013


Hi Fraser,

I just realized that my previous email was everything else than “self-explanatory” and on top of that - 
the workaround contained a typo…

Cheers,
Okko

Please find below the revised code including an example at the end:

library(oligo)
library(farms)

farms <-  function(object, background=FALSE, normalize=TRUE, subset=NULL, target="core", weight=.5, 
                     mu=0, weighted.mean=FALSE, laplacian=FALSE, robust=TRUE, 
                     correction=0, centering=c("median","mean"), spuriousCorrelation=0, verbose=TRUE )
{		
  target <- match.arg(target, c("core", "probeset"))
  featureInfo <- oligo:::stArrayPmInfo(object, target=target)
  theClass <- class(exprs(object))
  pmi <- featureInfo[["fid"]]
  pnVec <- as.character(featureInfo[["fsetid"]])
  if ("matrix" %in% theClass){
    pms <- exprs(object)[pmi,, drop=FALSE]
    dimnames(pms) <- NULL
    colnames(pms) <- sampleNames(object)
    theExprs <- basicFARMS(pms, pnVec, normalize, background, weight=weight, mu=mu, weighted.mean=weighted.mean, 
                           laplacian=laplacian, robust=robust, correction=correction, centering=centering, spuriousCorrelation=spuriousCorrelation)
    rm(pms)
  }
  else{
    stop("farms is currently not implemented for '", theClass, "' objects.")
  }
  
  
  colnames(theExprs$exprs) <- colnames(theExprs$se.exprs) <-  sampleNames(object)
  rownames(theExprs$exprs) <- rownames(theExprs$se.exprs) <- unique(pnVec)     
  
  out <- new("ExpressionSet",
             phenoData=phenoData(object),
             experimentData=experimentData(object),
             exprs=theExprs$exprs,
             se.exprs=theExprs$se.exprs,
             annotation=annotation(object),
             protocolData=protocolData(object))
  
  if (validObject(out)){
    return(out)
  }else{
    stop("Resulting object is invalid.")
  }
}




basicFARMS <- function(pmMat, pnVec, normalize, background, weight, mu, weighted.mean, 
                       laplacian, robust, correction, centering, spuriousCorrelation){
  pns <- unique(pnVec)
  nPn <- length(unique(pnVec))
  pnVec <- split(1:(length(pnVec)), pnVec)
  narray <- ncol(pmMat)
  if(background){
    pmMat <- oligo:::backgroundCorrect(pmMat)
  }
  if(normalize){
    pmMat <- oligo:::normalize(pmMat)
  }
  ini.mat <- matrix(0.5, nPn, narray )
  exprs.mat <- matrix(0, nPn, narray )
  message("Summarizing... ", appendLF=FALSE)
  for(i in 1:length(pnVec)){
    if(length(pnVec[[i]])<2){
      exprs.mat[i,] <- log2(pmMat[pnVec[[i]],])
    }
    else{
      res <- generateExprVal.method.farms(probes=pmMat[pnVec[[i]],], weight=weight, mu=mu,  cyc=30, tol=0.00001, 
                                          weighted.mean=weighted.mean, robust=robust, minNoise=minNoise, correction=correction,
                                          laplacian=laplacian, centering=centering, spuriousCorrelation=spuriousCorrelation)
      exprs.mat[i,] <- res$exprs
      ini.mat[i,] <- res$se.exprs
    }
  }
  message("OK")
  return(list(exprs=exprs.mat, se.exprs=ini.mat))
}



### example to run 
# library(oligo)
# library(farms)
# celFiles <- list.celfiles("/path/to/celfiles/", full.name=TRUE)
# raw <- read.celfiles(celFiles)
# eset <- farms(raw, normalize=TRUE)
####





-- 
Djork-Arné Clevert, PhD

Phone: +49 30 6883 5306
Fax: +49 30 6883 5307
Email: okko at clevert.de

Am 13.11.2013 um 14:32 schrieb Fraser Sim <fjsim at buffalo.edu>:

> Hi,
> 
> 
> 
> I am analyzing Affymetrix hugene2.0st arrays using the oligo package and had
> hoped to be able to use informative/non-informative probe calling (FARMS,
> Tolloen et 2007) as a means on filtering data. I had used this on older 3'
> biased arrays analyzed with affy but the farms package only works with an
> affyBatch object. Is this possible with another package or should I be
> considering a different approach to filter array data on the newer Affy
> arrays?
> 
> 
> 
> Cheers,
> 
> Fraser
> 
> 
> 	[[alternative HTML version deleted]]
> 
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor



More information about the Bioconductor mailing list