[BioC] trouble with gcrma background correction

[Zhijin (Jean) Wu]

HI, Sucheta,
   I haven't read your script but maybe it does the same as function 
"bg.adjust.gcrma" as defined in the gcrma package. Did you edit from that 
function and simply added the messages? Or is there any other change?
Jean

On Fri, 17 Nov 2006, Sucheta Tripathy wrote:

> Dear Group,
>
> I am stuck with this problem for quite some time. I have been searching
> for an answer to this, but could not quite get a hold of it. Any help in
> this will be greatly appreciated.
>
> I am using R version 2.3.1 and bioconductor version 1.9.
>
> First I computed affinity using the following command:
>
> affinity.info <- compute.affinities("Soybean", verbose=TRUE)
> save(affinity.info,file = "soybean.RData");
>
> Then I load it using:
>
> library(gcrma)
> load("soybean.RData")
> source("gcrmaBG.R")
>
> when I run the script gcrmaBG.R inside another script, I get the following
> error:
>
> ERROR MESSAGES:
> ===============
> Adjusting for optical effect.Done.
> In first if and type is
> Adjusting for non-specific binding
> .type is fullmodel
> Error in model.frame(formula, rownames, variables, varnames, extras,
> extranames,
>  :
>    invalid variable type for 'affinities'
> Execution halted
>
> ==========
> gcrmaBG.R
> ==========
>
> gcrmaBG <- function(object,
>                  type=c("fullmodel","affinities","mm","constant"),
>                  k=6*fast+0.5*(1-fast),stretch=1.15*fast+1*(1-fast),
>                  correction=1,rho=0.7,
>                  optical.correct=TRUE,verbose=TRUE,fast=TRUE){
>
>
>  type <- match.arg(type)
>
>
>  pmonly  <- (type=="affinities"|type=="constant")
>  needaff <- (type=="fullmodel"|type=="affinities")
>
>
>  if(optical.correct){
>    object <- bg.adjust.optical(object,verbose=verbose)
>  }
>
>
>  pms_bg <- gcrma.engine(pms=pm(object),
>                         mms=mm(object),
>                         pm.affinities=pm(affinity.info),
>                         mm.affinities=mm(affinity.info),
>                         type=type,k=k,
>                         stretch=stretch,
>                         correction=correction,rho=rho,
>                         verbose=verbose,fast=fast)
>
>  return(pms_bg)
> }
>
>
> bg.adjust.optical <- function(abatch,minimum=1,verbose=TRUE){
>  Index <- unlist(indexProbes(abatch,"both"))
>
>
>    if(verbose) cat("Adjusting for optical effect")
>
>        for(i in 1:length(abatch)){
>            if(verbose) cat(".")
>            exprs(abatch)[Index,i] <- exprs(abatch)[Index,i] -
>            min(exprs(abatch)[Index,i],na.rm=TRUE) + minimum
>        }
>
>    if(verbose) cat("Done.\n")
>
>    abatch
>
>    }
>
>
> # Defining function gcrma.engine
>
> gcrma.engine <- function(pms,mms,pm.affinities=NULL,mm.affinities=NULL,
>                         type=c("fullmodel","affinities","mm","constant"),
>                         k=6*fast+0.25*(1-fast),
>                         stretch=1.15*fast+1*(1-fast),correction=1,rho=0.7,
>                         verbose=TRUE,fast=TRUE){
>  type <- match.arg(type)
>
>  if(!is.null(pm.affinities)){
>    index.affinities <- which(!is.na(pm.affinities))
>    pm.affinities <- pm.affinities[index.affinities]
>
>    if(!is.null(mm.affinities)){
>      mm.affinities <- mm.affinities[index.affinities]
>    }
>
>  }
  >
>
>  if(type=="fullmodel" | type=="affinities"){
>
>  cat("In first if and type is \n")
>
>    set.seed(1)
>    Subset <- sample(1:length(pms[index.affinities,]),25000)
>    y <- log2(pms)[index.affinities,][Subset]
>    Subset <- (Subset-1)%%length(pms[index.affinities,])+1
>    x <- pm.affinities[Subset]
>    fit1 <- lm(y~x)
>  }
>
>
>  cat("Adjusting for non-specific binding\n")
>  for(i in 1:ncol(pms)){
>    if(verbose) cat(".")
>    if(type=="fullmodel"){
>
>    cat("type is fullmodel\n")  # ERROR HERE
>
>    pms[,i] <- bg.adjust.fullmodel(pms[,i],mms[,i],
>                                     pm.affinities,mm.affinities,
>                                     index.affinities,k=k,
>                                     rho=rho,fast=fast)
>    cat("bg.adjust.fullmodel is done\n")
>
>      pms[index.affinities,i] <- 2^(log2(pms[index.affinities,i])-
>                                    fit1$coef[2]*pm.affinities+mean(fit1$coef[2]
> *pm.affinities))
>    }
>
>    if(type=="affinities"){
>    cat("Type is affinities")
>      pms[,i] <- bg.adjust.affinities(pms[,i],mms[,i],
>                                      pm.affinities,mm.affinities,
>                                      index.affinities, k=k,
>                                      fast=fast)
>      pms[index.affinities,i] <- 2^(log2(pms[index.affinities,i])-
>                                    fit1$coef[2]*pm.affinities +
>                                    mean(fit1$coef[2]*pm.affinities))
>    }
>
>    if(type=="mm") pms[,i] <-
> bg.adjust.mm(pms[,i],correction*mms[,i],k=k,fast=f
> ast)
>
>        cat("type is mm") # added by ST.
>
>    if(type=="constant"){
>    cat("Type is constant")
>      pms[,i] <-
> bg.adjust.constant(pms[,i],k=k,Q=correction*mean(pms[,i]<mms[,i
> ]),fast=fast)
>    }
>    if(stretch!=1){
>      mu <- mean(log(pms[,i]))
>      pms[,i] <- exp(mu + stretch*(log(pms[,i])-mu))
>    }
>  }
>
>  if(verbose) cat("Done.\n")
>
>  return(pms)
> }
>
> Many thanks
>
> Sucheta
>
> -- 
> Sucheta Tripathy, Ph.D.
> Virginia Bioinformatics Institute Phase-I
> Washington street.
> Virginia Tech.
> Blacksburg,VA 24061-0447
> phone:(540)231-8138
> Fax:  (540) 231-2606
>
> web page: http://staff.vbi.vt.edu/sutripa
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor
>