[BioC] Trouble with gcrma background correction

[Herve Pages]

Hi Sucheta,

Sucheta Tripathy wrote:
> Dear Group,
> 
> (Sorry for the repost- The earlier one could not get through)
> 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.

We strongly discourage you of doing this. Bioconductor 1.9 has been
designed and tested for R 2.4.z. It is not supported for R 2.3.z users
so you are on your own...

> 
> 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


I suggest you upgrade to R 2.4.0 first and see if you can reproduce this
error. Also, don't forget to provide the output of sessionInfo() so, among
other things, people know with version of gcrma you are using.

One last advice: the long script you provided below doesn't help much
because:

  (1) some lines have been cut in the middle of a variable name
            if(type=="mm") pms[,i] <-
        bg.adjust.mm(pms[,i],correction*mms[,i],k=k,fast=f
        ast)
      so we can't just copy and paste it into a file and expect it
      to work without some extra editing work

  (2) it's too long

Please try to reproduce the error in a few lines of code (no more than 10)
and send only those 10 lines (with sessionInfo output).

Thanks,

H.


> 
> ==========
> 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
>