[BioC] Trouble with gcrma background correction

[Sucheta Tripathy]

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.

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
blog    : http://genomics-array.blogspot.com/