R/snm.R

In snm: Supervised Normalization of Microarrays

snm <-
	function(raw.dat, bio.var=NULL, adj.var=NULL, int.var=NULL,
                 weights=NULL, spline.dim = 4, num.iter = 10, lmer.max.iter=1000,
                 nbins=20, rm.adj=FALSE, verbose=TRUE, diagnose=TRUE)
{

  if(is.null(bio.var)) {
    warning("bio.var=NULL, so all probes will be treated as 'null' in the normalization.")
  }
  
  snm.obj <- make.snm.obj(Y=raw.dat, bio.var, adj.var, int.var, spline.dim,
      nbins, weights, diagnose, rm.adj)  

  if(!is.null(snm.obj$int.var)) {		  

    #Set up model fitting loop
    basisSplineFunction <- buildBasisFunction(snm.obj)
    
    #If bio.var=NULL, then only fit intensity-dependent effects
    if(is.null(snm.obj$bio.var)) { 
      snm.obj$pi0 <- 1
      snm.obj$nulls <- 1:snm.obj$n.probes
      obs.fit = list()
      xx0 = snm.obj$adj.var
      H0 = xx0 %*% solve(t(xx0) %*% xx0) %*% t(xx0)
      obs.fit$fit0 = t(H0 %*% t(snm.obj$dat))
      obs.fit$res1 = snm.obj$dat - obs.fit$fit0
      snm.obj <- fit.model(obs.fit, snm.obj, basisSplineFunction, lmer.max.iter)
      snm.obj$dat <- snm.obj$r.dat - snm.obj$array.fx
      snm.obj$pi0s <- NULL
    } else {
    #Otherwise perform model fitting loop
      pi0s <- rep(0,num.iter)
      obs.fit <- edge.fit(snm.obj, odp=FALSE)
      for (i in 1:num.iter) {
        if(verbose){cat("\r", "Iteration: ", i)}
        obs.stat <- edge.glr(obs.fit, df1=snm.obj$df.full, df0=snm.obj$df.null, norm.pval=TRUE)
      	snm.obj$pi0 <- edge.qvalue(obs.stat$pval)$pi0
      	snm.obj$nulls <- calculate.nulls(obs.stat$pval, snm.obj$pi0)
        ## These two steps should be included to follow the paper exactly
        ## snm.obj$dat <- snm.obj$r.dat
        ## obs.fit <- edge.fit(snm.obj, odp=FALSE)
        snm.obj <- fit.model(obs.fit, snm.obj, basisSplineFunction, lmer.max.iter)
        snm.obj$dat <- snm.obj$r.dat - snm.obj$array.fx
      	obs.fit <- edge.fit(snm.obj, odp=FALSE)
      	if(diagnose) {
          snm.diagnostic.plot(obs.fit,obs.stat$pval,snm.obj,iter=i)
      	}
      	pi0s[i] <- snm.obj$pi0
      }
      snm.obj$pi0s <- pi0s
    }

    #Extract subset of array effects (to reduce size of return object)
    #Set seed before running snm if this needs to be reproduced
    if(dim(snm.obj$dat)[1] > (nbins*100)) {
      th <- sample(dim(snm.obj$dat)[1], (nbins*100))
    } else {
      th <- 1:dim(snm.obj$dat)[1]
    } 
    M.ret = snm.obj$M[th,]
    array.fx.ret = snm.obj$array.fx[th,]
    for(i in 1:ncol(M.ret)) {
      oo <- order(M.ret[,i])
      M.ret[,i] = M.ret[oo,i]
      array.fx.ret[,i] = array.fx.ret[oo,i]
    }

    #Get final model fit significance
    if(is.null(bio.var)) { 
      obs.stat <- NULL
      snm.obj$pval <- NULL
      snm.obj$pi0 <- NULL
    } else {
      obs.stat <- edge.glr(obs.fit, df1=snm.obj$df.full, df0=snm.obj$df.null, norm.pval=TRUE)
      snm.obj$pval <- obs.stat$pval
      snm.obj$pi0 <- edge.qvalue(obs.stat$pval)$pi0
      if(snm.obj$rm.adj){
      	snm.obj <- remove.adjustment.vars(snm.obj)
      }
    }
  } else {
    #No intensity-dep effects, so simply remove adjustment variables
    obs.fit <- edge.fit(snm.obj, odp=FALSE)
    obs.stat <- edge.glr(obs.fit, df1=snm.obj$df.full, df0=snm.obj$df.null, norm.pval=TRUE)
    if(!snm.obj$rm.adj) {
      stop("int.var=NULL and rm.adj=FALSE, so there is nothing to do.")
    }
    snm.obj <- remove.adjustment.vars(snm.obj)
    snm.obj$pval <- obs.stat$pval; 
    snm.obj$pi0s <- edge.qvalue(obs.stat$pval)$pi0
    snm.obj$pi0 <- snm.obj$pi0s
    M.ret = NULL
    array.fx.ret = NULL
  }

  dimnames(snm.obj$adj.var) = snm.obj$dimnames.adj
  dimnames(snm.obj$bio.var) = snm.obj$dimnames.bio
	
  snm.ret <- list(norm.dat=snm.obj$dat, pval=snm.obj$pval, pi0=snm.obj$pi0, iter.pi0s=snm.obj$pi0s, nulls=snm.obj$nulls,
                  M=M.ret, array.fx=array.fx.ret, bio.var=snm.obj$bio.var, adj.var=snm.obj$adj.var, int.var=snm.obj$int.var,
                  df1=snm.obj$df.full, df0=snm.obj$df.null, raw.dat=raw.dat, rm.adj=rm.adj, call = match.call())

  class(snm.ret) = "snm"
  return(snm.ret)
}