@@ -188,7 +188,6 @@ snprma2 <- function(filenames, mixtureSampleSize=10^5, fitMixture=FALSE,

 processCEL <- function(i, filenames, fitMixture, A, B, SKW, SNR,

                        mixtureParams, eps, seed, mixtureSampleSize,

                        pkgname){

-

   obj1 <- loader("preprocStuff.rda", .crlmmPkgEnv, pkgname)

   obj2 <- loader("genotypeStuff.rda", .crlmmPkgEnv, pkgname)

   obj3 <- loader("mixtureStuff.rda", .crlmmPkgEnv, pkgname)

@@ -15,7 +15,7 @@ snprma <- function(filenames, mixtureSampleSize=10^5, fitMixture=FALSE,

     message(strwrap(msg))

     stop("Package ", pkgname, " could not be found.")

   }

-  

+

   if(verbose) message("Loading annotations and mixture model parameters.")

   loader("preprocStuff.rda", .crlmmPkgEnv, pkgname)

   loader("genotypeStuff.rda", .crlmmPkgEnv, pkgname)

@@ -30,7 +30,7 @@ snprma <- function(filenames, mixtureSampleSize=10^5, fitMixture=FALSE,

   SMEDIAN <- getVarInEnv("SMEDIAN")

   theKnots <- getVarInEnv("theKnots")

   gns <- getVarInEnv("gns")

-  

+

   ##We will read each cel file, summarize, and run EM one by one

   ##We will save parameters of EM to use later

   mixtureParams <- matrix(0, 4, length(filenames))

@@ -39,7 +39,7 @@ snprma <- function(filenames, mixtureSampleSize=10^5, fitMixture=FALSE,

 ##   mixtureParams <- initializeBigMatrix("crlmmMixt-", 4, length(filenames))

 ##   SNR <- initializeBigVector("crlmmSNR-", length(filenames), "numeric")

 ##   SKW <- initializeBigVector("crlmmSKW-", length(filenames), "numeric")

-  

+

   ## This is the sample for the fitting of splines

   ## BC: I like better the idea of the user passing the seed,

   ##     because this might intefere with other analyses

@@ -51,15 +51,15 @@ snprma <- function(filenames, mixtureSampleSize=10^5, fitMixture=FALSE,

   ##f is the correction. we save to avoid recomputing

   A <- matrix(as.integer(0), length(pnsa), length(filenames))

   B <- matrix(as.integer(0), length(pnsb), length(filenames))

-  

+

   if(verbose){

     message("Processing ", length(filenames), " files.")

     pb <- txtProgressBar(min=0, max=length(filenames), style=3)

   }

-  

+

   ##for skewness. no need to do everything

   idx2 <- sample(length(fid), 10^5)

-  

+

   ##We start looping throug cel files

   for(i in seq(along=filenames)){

     y <- as.matrix(read.celfile(filenames[i], intensity.means.only=TRUE)[["INTENSITY"]][["MEAN"]][fid])

@@ -73,10 +73,10 @@ snprma <- function(filenames, mixtureSampleSize=10^5, fitMixture=FALSE,

     if(fitMixture){

       S <- (log2(A[idx, i])+log2(B[idx, i]))/2 - SMEDIAN

       M <- log2(A[idx, i])-log2(B[idx, i])

-      

+

       ##we need to test the choice of eps.. it is not the max diff between funcs

       tmp <- fitAffySnpMixture56(S, M, theKnots, eps=eps)

-      

+

       mixtureParams[, i] <- tmp[["coef"]]

       SNR[i] <- tmp[["medF1"]]^2/(tmp[["sigma1"]]^2+tmp[["sigma2"]]^2)

     }

@@ -97,12 +97,12 @@ fitAffySnpMixture56 <- function(S, M, knots, probs=rep(1/3, 3), eps=.01, maxit=1

   mus <- append(quantile(M, c(1, 5)/6, names=FALSE), 0, 1)

   sigmas <- rep(mad(c(M[M<mus[1]]-mus[1], M[M>mus[3]]-mus[3])), 3)

   sigmas[2] <- sigmas[2]/2

- 

+

   weights <- apply(cbind(mus, sigmas), 1, function(p) dnorm(M, p[1], p[2]))

   previousF1 <- -Inf

   change <- eps+1

   it <- 0

- 

+

   if(verbose) message("Max change must be under ", eps, ".")

   matS <- stupidSplineBasis(S, knots)

   while (change > eps & it < maxit){

@@ -112,19 +112,19 @@ fitAffySnpMixture56 <- function(S, M, knots, probs=rep(1/3, 3), eps=.01, maxit=1

     LogLik <- rowSums(z)

     z <- sweep(z, 1, LogLik, "/")

     probs <- colMeans(z)

- 

+

     ## M

     fit1 <- crossprod(chol2inv(chol(crossprod(sweep(matS, 1, z[, 1], FUN="*"), matS))), crossprod(matS, z[, 1]*M))

- 

+

     fit2 <- sum(z[, 2]*M)/sum(z[, 2])

     F1 <- matS%*%fit1

     sigmas[c(1, 3)] <- sqrt(sum(z[, 1]*(M-F1)^2)/sum(z[, 1]))

     sigmas[2] <- sqrt(sum(z[, 2]*(M-fit2)^2)/sum(z[, 2]))

- 

+

     weights[, 1] <- dnorm(M, F1, sigmas[1])

     weights[, 2] <- dnorm(M, fit2, sigmas[2])

     weights[, 3] <- dnorm(M, -F1, sigmas[3])

-    

+

     change <- max(abs(F1-previousF1))

     previousF1 <- F1

     if(verbose) message("Iter ", it, ": ", change, ".")

@@ -140,7 +140,7 @@ snprma2 <- function(filenames, mixtureSampleSize=10^5, fitMixture=FALSE,

     cdfName <- read.celfile.header(filenames[1])[["cdfName"]]

   pkgname <- getCrlmmAnnotationName(cdfName)

   stopifnot(require(pkgname, character.only=TRUE, quietly=!verbose))

-  

+

   if(verbose) message("Loading annotations and mixture model parameters.")

   obj <- loader("preprocStuff.rda", .crlmmPkgEnv, pkgname)

   pnsa <- getVarInEnv("pnsa")

@@ -148,14 +148,14 @@ snprma2 <- function(filenames, mixtureSampleSize=10^5, fitMixture=FALSE,

   gns <- getVarInEnv("gns")

   rm(list=obj, envir=.crlmmPkgEnv)

   rm(obj)

-  

+

   ##We will read each cel file, summarize, and run EM one by one

   ##We will save parameters of EM to use later

   if(verbose) message("Initializing objects.")

   mixtureParams <- initializeBigMatrix("crlmmMixt-", 4, length(filenames), "double")

   SNR <- initializeBigVector("crlmmSNR-", length(filenames), "double")

   SKW <- initializeBigVector("crlmmSKW-", length(filenames), "double")

-  

+

   ## This is the sample for the fitting of splines

   ## BC: I like better the idea of the user passing the seed,

   ##     because this might intefere with other analyses

@@ -180,7 +180,7 @@ snprma2 <- function(filenames, mixtureSampleSize=10^5, fitMixture=FALSE,

   close(SKW)

   close(A)

   close(B)

-  

+

   list(A=A, B=B, sns=sns, gns=gns, SNR=SNR, SKW=SKW, mixtureParams=mixtureParams, cdfName=cdfName)

 }

@@ -188,7 +188,7 @@ snprma2 <- function(filenames, mixtureSampleSize=10^5, fitMixture=FALSE,

 processCEL <- function(i, filenames, fitMixture, A, B, SKW, SNR,

                        mixtureParams, eps, seed, mixtureSampleSize,

                        pkgname){

-  

+

   obj1 <- loader("preprocStuff.rda", .crlmmPkgEnv, pkgname)

   obj2 <- loader("genotypeStuff.rda", .crlmmPkgEnv, pkgname)

   obj3 <- loader("mixtureStuff.rda", .crlmmPkgEnv, pkgname)

@@ -226,7 +226,7 @@ processCEL <- function(i, filenames, fitMixture, A, B, SKW, SNR,

     A[, k] <- intMedianSummaries(y[aIndex, 1, drop=FALSE], pnsa)

     B[, k] <- intMedianSummaries(y[bIndex, 1, drop=FALSE], pnsb)

     rm(y)

-    

+

     if(fitMixture){

       S <- (log2(A[idx,k])+log2(B[idx, k]))/2 - SMEDIAN

       M <- log2(A[idx, k])-log2(B[idx, k])

@@ -142,10 +142,12 @@ snprma2 <- function(filenames, mixtureSampleSize=10^5, fitMixture=FALSE,

   stopifnot(require(pkgname, character.only=TRUE, quietly=!verbose))

   if(verbose) message("Loading annotations and mixture model parameters.")

-  loader("preprocStuff.rda", .crlmmPkgEnv, pkgname)

+  obj <- loader("preprocStuff.rda", .crlmmPkgEnv, pkgname)

   pnsa <- getVarInEnv("pnsa")

   pnsb <- getVarInEnv("pnsb")

   gns <- getVarInEnv("gns")

+  rm(list=obj, envir=.crlmmPkgEnv)

+  rm(obj)

   ##We will read each cel file, summarize, and run EM one by one

   ##We will save parameters of EM to use later

@@ -173,6 +175,11 @@ snprma2 <- function(filenames, mixtureSampleSize=10^5, fitMixture=FALSE,

            mixtureParams=mixtureParams, eps=eps, seed=seed,

            mixtureSampleSize=mixtureSampleSize, pkgname=pkgname,

            neededPkgs=c("crlmm", pkgname))

+  close(mixtureParams)

+  close(SNR)

+  close(SKW)

+  close(A)

+  close(B)

   list(A=A, B=B, sns=sns, gns=gns, SNR=SNR, SKW=SKW, mixtureParams=mixtureParams, cdfName=cdfName)

 }

@@ -182,9 +189,9 @@ processCEL <- function(i, filenames, fitMixture, A, B, SKW, SNR,

                        mixtureParams, eps, seed, mixtureSampleSize,

                        pkgname){

-  loader("preprocStuff.rda", .crlmmPkgEnv, pkgname)

-  loader("genotypeStuff.rda", .crlmmPkgEnv, pkgname)

-  loader("mixtureStuff.rda", .crlmmPkgEnv, pkgname)

+  obj1 <- loader("preprocStuff.rda", .crlmmPkgEnv, pkgname)

+  obj2 <- loader("genotypeStuff.rda", .crlmmPkgEnv, pkgname)

+  obj3 <- loader("mixtureStuff.rda", .crlmmPkgEnv, pkgname)

   autosomeIndex <- getVarInEnv("autosomeIndex")

   pnsa <- getVarInEnv("pnsa")

   pnsb <- getVarInEnv("pnsb")

@@ -195,6 +202,8 @@ processCEL <- function(i, filenames, fitMixture, A, B, SKW, SNR,

   SMEDIAN <- getVarInEnv("SMEDIAN")

   theKnots <- getVarInEnv("theKnots")

   gns <- getVarInEnv("gns")

+  rm(list=c(obj1, obj2, obj3), envir=.crlmmPkgEnv)

+  rm(obj1, obj2, obj3)

   ## for mixture

   set.seed(seed)

@@ -222,8 +231,10 @@ processCEL <- function(i, filenames, fitMixture, A, B, SKW, SNR,

       S <- (log2(A[idx,k])+log2(B[idx, k]))/2 - SMEDIAN

       M <- log2(A[idx, k])-log2(B[idx, k])

       tmp <- fitAffySnpMixture56(S, M, theKnots, eps=eps)

+      rm(S, M)

       mixtureParams[, k] <- tmp[["coef"]]

       SNR[k] <- tmp[["medF1"]]^2/(tmp[["sigma1"]]^2+tmp[["sigma2"]]^2)

+      rm(tmp)

     } else {

       mixtureParams[, k] <- NA

       SNR[k] <- NA

@@ -234,5 +245,7 @@ processCEL <- function(i, filenames, fitMixture, A, B, SKW, SNR,

   close(SKW)

   close(mixtureParams)

   close(SNR)

+  rm(list=ls())

+  gc(verbose=FALSE)

   TRUE

 }

@@ -1,17 +1,21 @@

-snprma <- function(filenames, mixtureSampleSize=10^5, fitMixture=FALSE, eps=0.1, verbose=TRUE, seed=1, cdfName, sns){

+snprma <- function(filenames, mixtureSampleSize=10^5, fitMixture=FALSE,

+                   eps=0.1, verbose=TRUE, seed=1, cdfName, sns){

   if (missing(sns)) sns <- basename(filenames)

-  ##ADD CHECK TO SEE IF LOADED

   if (missing(cdfName))

     cdfName <- read.celfile.header(filenames[1])$cdfName

-  ##  stuffDir <- changeToCrlmmAnnotationName(cdfName)

   pkgname <- getCrlmmAnnotationName(cdfName)

+

   if(!require(pkgname, character.only=TRUE, quietly=!verbose)){

-    suggCall <- paste("library(", pkgname, ", lib.loc='/Altern/Lib/Loc')", sep="")

-    msg <- paste("If", pkgname, "is installed on an alternative location, please load it manually by using", suggCall)

+    suggCall <- paste("library(", pkgname,

+                      ", lib.loc='/Altern/Lib/Loc')",

+                      sep="")

+    msg <- paste("If", pkgname,

+                 "is installed on an alternative location,",

+                 "please load it manually by using", suggCall)

     message(strwrap(msg))

     stop("Package ", pkgname, " could not be found.")

-    rm(suggCall, msg)

   }

+  

   if(verbose) message("Loading annotations and mixture model parameters.")

   loader("preprocStuff.rda", .crlmmPkgEnv, pkgname)

   loader("genotypeStuff.rda", .crlmmPkgEnv, pkgname)

@@ -32,6 +36,9 @@ snprma <- function(filenames, mixtureSampleSize=10^5, fitMixture=FALSE, eps=0.1,

   mixtureParams <- matrix(0, 4, length(filenames))

   SNR <- vector("numeric", length(filenames))

   SKW <- vector("numeric", length(filenames))

+##   mixtureParams <- initializeBigMatrix("crlmmMixt-", 4, length(filenames))

+##   SNR <- initializeBigVector("crlmmSNR-", length(filenames), "numeric")

+##   SKW <- initializeBigVector("crlmmSKW-", length(filenames), "numeric")

   ## This is the sample for the fitting of splines

   ## BC: I like better the idea of the user passing the seed,

@@ -47,10 +54,13 @@ snprma <- function(filenames, mixtureSampleSize=10^5, fitMixture=FALSE, eps=0.1,

   if(verbose){

     message("Processing ", length(filenames), " files.")

-    if (getRversion() > '2.7.0') pb <- txtProgressBar(min=0, max=length(filenames), style=3)

+    pb <- txtProgressBar(min=0, max=length(filenames), style=3)

   }

+  

+  ##for skewness. no need to do everything

+  idx2 <- sample(length(fid), 10^5)

+  

   ##We start looping throug cel files

-  idx2 <- sample(length(fid), 10^5) ##for skewness. no need to do everything

   for(i in seq(along=filenames)){

     y <- as.matrix(read.celfile(filenames[i], intensity.means.only=TRUE)[["INTENSITY"]][["MEAN"]][fid])

     x <- log2(y[idx2])

@@ -70,13 +80,9 @@ snprma <- function(filenames, mixtureSampleSize=10^5, fitMixture=FALSE, eps=0.1,

       mixtureParams[, i] <- tmp[["coef"]]

       SNR[i] <- tmp[["medF1"]]^2/(tmp[["sigma1"]]^2+tmp[["sigma2"]]^2)

     }

-    if (verbose)

-      if (getRversion() > '2.7.0') setTxtProgressBar(pb, i)

-      else cat(".")

+    if (verbose) setTxtProgressBar(pb, i)

   }

-  if (verbose)

-    if (getRversion() > '2.7.0') close(pb)

-    else cat("\n")

+  if (verbose) close(pb)

   if (!fitMixture) SNR <- mixtureParams <- NA

   ## gns comes from preprocStuff.rda

   list(A=A, B=B, sns=sns, gns=gns, SNR=SNR, SKW=SKW, mixtureParams=mixtureParams, cdfName=cdfName)

@@ -127,5 +133,106 @@ fitAffySnpMixture56 <- function(S, M, knots, probs=rep(1/3, 3), eps=.01, maxit=1

  return(list(coef= -fit1, medF1=medF1, sigma1=sigmas[1], sigma2=sigmas[2]))

 }

+snprma2 <- function(filenames, mixtureSampleSize=10^5, fitMixture=FALSE,

+                    eps=0.1, verbose=TRUE, seed=1, cdfName, sns){

+  if (missing(sns)) sns <- basename(filenames)

+  if (missing(cdfName))

+    cdfName <- read.celfile.header(filenames[1])[["cdfName"]]

+  pkgname <- getCrlmmAnnotationName(cdfName)

+  stopifnot(require(pkgname, character.only=TRUE, quietly=!verbose))

+  

+  if(verbose) message("Loading annotations and mixture model parameters.")

+  loader("preprocStuff.rda", .crlmmPkgEnv, pkgname)

+  pnsa <- getVarInEnv("pnsa")

+  pnsb <- getVarInEnv("pnsb")

+  gns <- getVarInEnv("gns")

+  

+  ##We will read each cel file, summarize, and run EM one by one

+  ##We will save parameters of EM to use later

+  if(verbose) message("Initializing objects.")

+  mixtureParams <- initializeBigMatrix("crlmmMixt-", 4, length(filenames), "double")

+  SNR <- initializeBigVector("crlmmSNR-", length(filenames), "double")

+  SKW <- initializeBigVector("crlmmSKW-", length(filenames), "double")

+  

+  ## This is the sample for the fitting of splines

+  ## BC: I like better the idea of the user passing the seed,

+  ##     because this might intefere with other analyses

+  ##     (like what happened to GCRMA)

+  ##S will hold (A+B)/2 and M will hold A-B

+  ##NOTE: We actually dont need to save S. Only for pics etc...

+  ##f is the correction. we save to avoid recomputing

+  A <- initializeBigMatrix("crlmmA-", length(pnsa), length(filenames), "integer")

+  B <- initializeBigMatrix("crlmmB-", length(pnsb), length(filenames), "integer")

+

+  sampleBatches <- splitIndicesByNode(seq(along=filenames))

+

+  if(verbose) message("Processing ", length(filenames), " files.")

+

+  ocLapply(sampleBatches, processCEL, filenames=filenames,

+           fitMixture=fitMixture, A=A, B=B, SKW=SKW, SNR=SNR,

+           mixtureParams=mixtureParams, eps=eps, seed=seed,

+           mixtureSampleSize=mixtureSampleSize, pkgname=pkgname,

+           neededPkgs=c("crlmm", pkgname))

+  

+  list(A=A, B=B, sns=sns, gns=gns, SNR=SNR, SKW=SKW, mixtureParams=mixtureParams, cdfName=cdfName)

+}

+

+processCEL <- function(i, filenames, fitMixture, A, B, SKW, SNR,

+                       mixtureParams, eps, seed, mixtureSampleSize,

+                       pkgname){

+  

+  loader("preprocStuff.rda", .crlmmPkgEnv, pkgname)

+  loader("genotypeStuff.rda", .crlmmPkgEnv, pkgname)

+  loader("mixtureStuff.rda", .crlmmPkgEnv, pkgname)

+  autosomeIndex <- getVarInEnv("autosomeIndex")

+  pnsa <- getVarInEnv("pnsa")

+  pnsb <- getVarInEnv("pnsb")

+  fid <- getVarInEnv("fid")

+  reference <- getVarInEnv("reference")

+  aIndex <- getVarInEnv("aIndex")

+  bIndex <- getVarInEnv("bIndex")

+  SMEDIAN <- getVarInEnv("SMEDIAN")

+  theKnots <- getVarInEnv("theKnots")

+  gns <- getVarInEnv("gns")

+

+  ## for mixture

+  set.seed(seed)

+  idx <- sort(sample(autosomeIndex, mixtureSampleSize))

+  ##for skewness. no need to do everything

+  idx2 <- sample(length(fid), 10^5)

+  open(A)

+  open(B)

+  open(SKW)

+  open(mixtureParams)

+  open(SNR)

+

+  for (k in i){

+    y <- as.matrix(read.celfile(filenames[k], intensity.means.only=TRUE)[["INTENSITY"]][["MEAN"]][fid])

+    x <- log2(y[idx2])

+    SKW[k] <- mean((x-mean(x))^3)/(sd(x)^3)

+    rm(x)

+    y <- normalize.quantiles.use.target(y, target=reference)

+    A[, k] <- intMedianSummaries(y[aIndex, 1, drop=FALSE], pnsa)

+    B[, k] <- intMedianSummaries(y[bIndex, 1, drop=FALSE], pnsb)

+    rm(y)

+    

+    if(fitMixture){

+      S <- (log2(A[idx,k])+log2(B[idx, k]))/2 - SMEDIAN

+      M <- log2(A[idx, k])-log2(B[idx, k])

+      tmp <- fitAffySnpMixture56(S, M, theKnots, eps=eps)

+      mixtureParams[, k] <- tmp[["coef"]]

+      SNR[k] <- tmp[["medF1"]]^2/(tmp[["sigma1"]]^2+tmp[["sigma2"]]^2)

+    } else {

+      mixtureParams[, k] <- NA

+      SNR[k] <- NA

+    }

+  }

+  close(A)

+  close(B)

+  close(SKW)

+  close(mixtureParams)

+  close(SNR)

+  TRUE

+}

@@ -127,3 +127,5 @@ fitAffySnpMixture56 <- function(S, M, knots, probs=rep(1/3, 3), eps=.01, maxit=1

  return(list(coef= -fit1, medF1=medF1, sigma1=sigmas[1], sigma2=sigmas[2]))

 }

+

+

@@ -42,11 +42,8 @@ snprma <- function(filenames, mixtureSampleSize=10^5, fitMixture=FALSE, eps=0.1,

   ##S will hold (A+B)/2 and M will hold A-B

   ##NOTE: We actually dont need to save S. Only for pics etc...

   ##f is the correction. we save to avoid recomputing

-  ##**RS**

-  ##A <- matrix(as.integer(0), length(pnsa), length(filenames))

-  ##B <- matrix(as.integer(0), length(pnsb), length(filenames))

-  A <- initializeBigMatrix(name="tmpA", length(pnsa), length(filenames))

-  B <- initializeBigMatrix(name="tmpB", length(pnsa), length(filenames))

+  A <- matrix(as.integer(0), length(pnsa), length(filenames))

+  B <- matrix(as.integer(0), length(pnsb), length(filenames))

   if(verbose){

     message("Processing ", length(filenames), " files.")

@@ -42,8 +42,11 @@ snprma <- function(filenames, mixtureSampleSize=10^5, fitMixture=FALSE, eps=0.1,

   ##S will hold (A+B)/2 and M will hold A-B

   ##NOTE: We actually dont need to save S. Only for pics etc...

   ##f is the correction. we save to avoid recomputing

-  A <- matrix(as.integer(0), length(pnsa), length(filenames))

-  B <- matrix(as.integer(0), length(pnsb), length(filenames))

+  ##**RS**

+  ##A <- matrix(as.integer(0), length(pnsa), length(filenames))

+  ##B <- matrix(as.integer(0), length(pnsb), length(filenames))

+  A <- initializeBigMatrix(name="tmpA", length(pnsa), length(filenames))

+  B <- initializeBigMatrix(name="tmpB", length(pnsa), length(filenames))

   if(verbose){

     message("Processing ", length(filenames), " files.")

@@ -1,129 +1,85 @@

-setMethod("snprma", "AffymetrixAlleleSet", function(object, filenames){

-	snprma.AffymetrixAlleleSet(object, filenames)

-})

-setMethod("snprma", "NChannelSet", function(object, ops, ...){

-	object <- preprocessInfinium2(object, ops)

-})

-

-snprma.AffymetrixAlleleSet <- function(object, filenames){

-##		   filenames,

-##		   mixtureSampleSize=10^5,

-##		   fitMixture=FALSE,

-##		   eps=0.1,

-##		   verbose=TRUE,

-##		   seed=1,

-##		   cdfName,

-##		   AFile="A.rda",

-##		   BFile="B.rda"){

-	ops <- crlmmOptions(object)

-	rmaOps <- ops$snprmaOpts

-	mixtureSampleSize <- rmaOps$mixtureSampleSize

-	fitMixture <- rmaOps$fitMixture

-	eps <- rmaOps$eps

-	verbose <- ops$verbose

-	seed <- ops$seed

-	cdfName <- annotation(object)

-	A <- A(object)

-	B <- B(object)

-	sns <- basename(filenames)

-	##ADD CHECK TO SEE IF LOADED

-	if (missing(cdfName))

-		cdfName <- read.celfile.header(filenames[1])$cdfName

-	##  stuffDir <- changeToCrlmmAnnotationName(cdfName)

-	pkgname <- getCrlmmAnnotationName(cdfName)

-	if(!require(pkgname, character.only=TRUE, quietly=!verbose)){

-		suggCall <- paste("library(", pkgname, ", lib.loc='/Altern/Lib/Loc')", sep="")

-		msg <- paste("If", pkgname, "is installed on an alternative location, please load it manually by using", suggCall)

-		message(strwrap(msg))

-		stop("Package ", pkgname, " could not be found.")

-		rm(suggCall, msg)

-	}

-	if(verbose) message("Loading annotations and mixture model parameters.")

-	loader("preprocStuff.rda", .crlmmPkgEnv, pkgname)

-	loader("genotypeStuff.rda", .crlmmPkgEnv, pkgname)

-	loader("mixtureStuff.rda", .crlmmPkgEnv, pkgname)

-	autosomeIndex <- getVarInEnv("autosomeIndex")

-	pnsa <- getVarInEnv("pnsa")

-	pnsb <- getVarInEnv("pnsb")

-	fid <- getVarInEnv("fid")

-	reference <- getVarInEnv("reference")

-	aIndex <- getVarInEnv("aIndex")

-	bIndex <- getVarInEnv("bIndex")

-	SMEDIAN <- getVarInEnv("SMEDIAN")

-	theKnots <- getVarInEnv("theKnots")

-	gns <- getVarInEnv("gns")

-	snpIndex <- which(isSnp(object)==1)

-	stopifnot(identical(featureNames(object)[snpIndex], gns))

+snprma <- function(filenames, mixtureSampleSize=10^5, fitMixture=FALSE, eps=0.1, verbose=TRUE, seed=1, cdfName, sns){

+  if (missing(sns)) sns <- basename(filenames)

+  ##ADD CHECK TO SEE IF LOADED

+  if (missing(cdfName))

+    cdfName <- read.celfile.header(filenames[1])$cdfName

+  ##  stuffDir <- changeToCrlmmAnnotationName(cdfName)

+  pkgname <- getCrlmmAnnotationName(cdfName)

+  if(!require(pkgname, character.only=TRUE, quietly=!verbose)){

+    suggCall <- paste("library(", pkgname, ", lib.loc='/Altern/Lib/Loc')", sep="")

+    msg <- paste("If", pkgname, "is installed on an alternative location, please load it manually by using", suggCall)

+    message(strwrap(msg))

+    stop("Package ", pkgname, " could not be found.")

+    rm(suggCall, msg)

+  }

+  if(verbose) message("Loading annotations and mixture model parameters.")

+  loader("preprocStuff.rda", .crlmmPkgEnv, pkgname)

+  loader("genotypeStuff.rda", .crlmmPkgEnv, pkgname)

+  loader("mixtureStuff.rda", .crlmmPkgEnv, pkgname)

+  autosomeIndex <- getVarInEnv("autosomeIndex")

+  pnsa <- getVarInEnv("pnsa")

+  pnsb <- getVarInEnv("pnsb")

+  fid <- getVarInEnv("fid")

+  reference <- getVarInEnv("reference")

+  aIndex <- getVarInEnv("aIndex")

+  bIndex <- getVarInEnv("bIndex")

+  SMEDIAN <- getVarInEnv("SMEDIAN")

+  theKnots <- getVarInEnv("theKnots")

+  gns <- getVarInEnv("gns")

+  

+  ##We will read each cel file, summarize, and run EM one by one

+  ##We will save parameters of EM to use later

+  mixtureParams <- matrix(0, 4, length(filenames))

+  SNR <- vector("numeric", length(filenames))

+  SKW <- vector("numeric", length(filenames))

-	##We will read each cel file, summarize, and run EM one by one

-	##We will save parameters of EM to use later

-	mixtureParams <- matrix(0, 4, length(filenames))

-	colnames(mixtureParams) <- basename(filenames)

-	SNR <- vector("numeric", length(filenames))

-	SKW <- vector("numeric", length(filenames))

+  ## This is the sample for the fitting of splines

+  ## BC: I like better the idea of the user passing the seed,

+  ##     because this might intefere with other analyses

+  ##     (like what happened to GCRMA)

+  set.seed(seed)

+  idx <- sort(sample(autosomeIndex, mixtureSampleSize))

+  ##S will hold (A+B)/2 and M will hold A-B

+  ##NOTE: We actually dont need to save S. Only for pics etc...

+  ##f is the correction. we save to avoid recomputing

+  A <- matrix(as.integer(0), length(pnsa), length(filenames))

+  B <- matrix(as.integer(0), length(pnsb), length(filenames))

-	## This is the sample for the fitting of splines

-	## BC: I like better the idea of the user passing the seed,

-	##     because this might intefere with other analyses

-	##     (like what happened to GCRMA)

-	set.seed(seed)

-	idx <- sort(sample(autosomeIndex, mixtureSampleSize))

-	##S will hold (A+B)/2 and M will hold A-B

-	##NOTE: We actually dont need to save S. Only for pics etc...

-	##f is the correction. we save to avoid recomputing

-	if(verbose){

-		message("Processing ", length(filenames), " files.")

-		if (getRversion() > '2.7.0') pb <- txtProgressBar(min=0, max=length(filenames), style=3)

-	}

-	##We start looping through cel files

-	idx2 <- sample(length(fid), 10^5) ##for skewness. no need to do everything

-	message("Quantile normalizing to HapMap target distribution")

-	for(i in seq(along=filenames)){

-		y <- as.matrix(read.celfile(filenames[i], intensity.means.only=TRUE)[["INTENSITY"]][["MEAN"]][fid])

-		x <- log2(y[idx2])

-		SKW[i] <- mean((x-mean(x))^3)/(sd(x)^3)

-		rm(x)

-		y <- normalize.quantiles.use.target(y, target=reference)

-		A[1:length(pnsa), i] <- intMedianSummaries(y[aIndex, 1, drop=FALSE], pnsa)

-		B[1:length(pnsa), i] <- intMedianSummaries(y[bIndex, 1, drop=FALSE], pnsb)

-		##Now to fit the EM

-		if(fitMixture){

-			S <- (log2(A[idx, i])+log2(B[idx, i]))/2 - SMEDIAN

-			M <- log2(A[idx, i])-log2(B[idx, i])

-			

-			##we need to test the choice of eps.. it is not the max diff between funcs

-			tmp <- fitAffySnpMixture56(S, M, theKnots, eps=eps)

-			

-			mixtureParams[, i] <- tmp[["coef"]]

-			SNR[i] <- tmp[["medF1"]]^2/(tmp[["sigma1"]]^2+tmp[["sigma2"]]^2)

-		}

-		if (verbose)

-			if (getRversion() > '2.7.0') setTxtProgressBar(pb, i)

-			else cat(".")

-	}

-##	A(object) <- A

-##	B(object) <- B

-	object@assayData[["alleleA"]] <- A

-	object@assayData[["alleleB"]] <- B

-	##closeObject(object, A)

-	##closeObject(object, B)

-	##if(isPackageLoaded("ff")) { close(A); close(B)} else {save(A, file=AFile); save(B, file=BFile)}

-	if (verbose)

-		if (getRversion() > '2.7.0') close(pb) else cat("\n")

-	##else cat("\n")

-	if (!fitMixture) SNR <- mixtureParams <- NA

-	sampleStats <- data.frame(SKW=SKW,

-				  SNR=SNR)

-##				  batch=ops$cnOpts$batch)

-	crlmmOpts <- crlmmOptions(object)$crlmmOpts

-	crlmmOpts$mixtureParams <- mixtureParams

-	crlmmOptions(object)$crlmmOpts <- crlmmOpts

-	pD <- new("AnnotatedDataFrame",

-		  data=sampleStats,

-		  varMetadata=data.frame(labelDescription=colnames(sampleStats)))

-	sampleNames(pD) <- sampleNames(object)

-	phenoData(object) <- pD

-	return(object)

+  if(verbose){

+    message("Processing ", length(filenames), " files.")

+    if (getRversion() > '2.7.0') pb <- txtProgressBar(min=0, max=length(filenames), style=3)

+  }

+  ##We start looping throug cel files

+  idx2 <- sample(length(fid), 10^5) ##for skewness. no need to do everything

+  for(i in seq(along=filenames)){

+    y <- as.matrix(read.celfile(filenames[i], intensity.means.only=TRUE)[["INTENSITY"]][["MEAN"]][fid])

+    x <- log2(y[idx2])

+    SKW[i] <- mean((x-mean(x))^3)/(sd(x)^3)

+    rm(x)

+    y <- normalize.quantiles.use.target(y, target=reference)

+    A[, i] <- intMedianSummaries(y[aIndex, 1, drop=FALSE], pnsa)

+    B[, i] <- intMedianSummaries(y[bIndex, 1, drop=FALSE], pnsb)

+    ##Now to fit the EM

+    if(fitMixture){

+      S <- (log2(A[idx, i])+log2(B[idx, i]))/2 - SMEDIAN

+      M <- log2(A[idx, i])-log2(B[idx, i])

+      

+      ##we need to test the choice of eps.. it is not the max diff between funcs

+      tmp <- fitAffySnpMixture56(S, M, theKnots, eps=eps)

+      

+      mixtureParams[, i] <- tmp[["coef"]]

+      SNR[i] <- tmp[["medF1"]]^2/(tmp[["sigma1"]]^2+tmp[["sigma2"]]^2)

+    }

+    if (verbose)

+      if (getRversion() > '2.7.0') setTxtProgressBar(pb, i)

+      else cat(".")

+  }

+  if (verbose)

+    if (getRversion() > '2.7.0') close(pb)

+    else cat("\n")

+  if (!fitMixture) SNR <- mixtureParams <- NA

+  ## gns comes from preprocStuff.rda

+  list(A=A, B=B, sns=sns, gns=gns, SNR=SNR, SKW=SKW, mixtureParams=mixtureParams, cdfName=cdfName)

 }

 fitAffySnpMixture56 <- function(S, M, knots, probs=rep(1/3, 3), eps=.01, maxit=10, verbose=FALSE){

@@ -1,7 +1,31 @@

-snprma <- function(filenames, mixtureSampleSize=10^5,

-		   fitMixture=FALSE, eps=0.1, verbose=TRUE,

-		   seed=1, cdfName, sns, AFile="A.rda", BFile="B.rda"){

-	if (missing(sns)) sns <- basename(filenames)

+setMethod("snprma", "AffymetrixAlleleSet", function(object, filenames){

+	snprma.AffymetrixAlleleSet(object, filenames)

+})

+setMethod("snprma", "NChannelSet", function(object, ops, ...){

+	object <- preprocessInfinium2(object, ops)

+})

+

+snprma.AffymetrixAlleleSet <- function(object, filenames){

+##		   filenames,

+##		   mixtureSampleSize=10^5,

+##		   fitMixture=FALSE,

+##		   eps=0.1,

+##		   verbose=TRUE,

+##		   seed=1,

+##		   cdfName,

+##		   AFile="A.rda",

+##		   BFile="B.rda"){

+	ops <- crlmmOptions(object)

+	rmaOps <- ops$snprmaOpts

+	mixtureSampleSize <- rmaOps$mixtureSampleSize

+	fitMixture <- rmaOps$fitMixture

+	eps <- rmaOps$eps

+	verbose <- ops$verbose

+	seed <- ops$seed

+	cdfName <- annotation(object)

+	A <- A(object)

+	B <- B(object)

+	sns <- basename(filenames)

 	##ADD CHECK TO SEE IF LOADED

 	if (missing(cdfName))

 		cdfName <- read.celfile.header(filenames[1])$cdfName

@@ -28,10 +52,13 @@ snprma <- function(filenames, mixtureSampleSize=10^5,

 	SMEDIAN <- getVarInEnv("SMEDIAN")

 	theKnots <- getVarInEnv("theKnots")

 	gns <- getVarInEnv("gns")

+	snpIndex <- which(isSnp(object)==1)

+	stopifnot(identical(featureNames(object)[snpIndex], gns))

 	##We will read each cel file, summarize, and run EM one by one

 	##We will save parameters of EM to use later

 	mixtureParams <- matrix(0, 4, length(filenames))

+	colnames(mixtureParams) <- basename(filenames)

 	SNR <- vector("numeric", length(filenames))

 	SKW <- vector("numeric", length(filenames))

@@ -44,23 +71,6 @@ snprma <- function(filenames, mixtureSampleSize=10^5,

 	##S will hold (A+B)/2 and M will hold A-B

 	##NOTE: We actually dont need to save S. Only for pics etc...

 	##f is the correction. we save to avoid recomputing

-  

-	path <- system.file("extdata", package=paste(cdfName, "Crlmm", sep=""))

-	load(file.path(path, "cnProbes.rda"))

-	cnProbes <- get("cnProbes")

-	nr <- length(pnsa) + nrow(cnProbes)

-	nc <- length(filenames)

-	dns <- list(c(gns, rownames(cnProbes)), sns)

-	if(isPackageLoaded("ff")){

-		message("ff package is loaded.  Creating ff objects for storing normalized and summarized intensities")

-		A <- initializeBigMatrix(nr, nc)

-		B <- initializeBigMatrix(nr, nc)

-		open(A)

-		open(B)

-	}  else {

-		A <- matrix(0L, nr, nc, dimnames=dns)

-		B <- matrix(0L, nr, nc, dimnames=dns)

-	}

 	if(verbose){

 		message("Processing ", length(filenames), " files.")

 		if (getRversion() > '2.7.0') pb <- txtProgressBar(min=0, max=length(filenames), style=3)

@@ -91,15 +101,29 @@ snprma <- function(filenames, mixtureSampleSize=10^5,

 			if (getRversion() > '2.7.0') setTxtProgressBar(pb, i)

 			else cat(".")

 	}

-	if(isPackageLoaded("ff")) { close(A); close(B)}

+##	A(object) <- A

+##	B(object) <- B

+	object@assayData[["alleleA"]] <- A

+	object@assayData[["alleleB"]] <- B

+	##closeObject(object, A)

+	##closeObject(object, B)

+	##if(isPackageLoaded("ff")) { close(A); close(B)} else {save(A, file=AFile); save(B, file=BFile)}

 	if (verbose)

-		if (getRversion() > '2.7.0') close(pb)

-		else cat("\n")

+		if (getRversion() > '2.7.0') close(pb) else cat("\n")

+	##else cat("\n")

 	if (!fitMixture) SNR <- mixtureParams <- NA

-	save(A, file=AFile)

-	save(B, file=BFile)

-	return(list(sns=sns, gns=gns, SNR=SNR, SKW=SKW, mixtureParams=mixtureParams, cdfName=cdfName,

-		    cnnames=rownames(cnProbes)))

+	sampleStats <- data.frame(SKW=SKW,

+				  SNR=SNR)

+##				  batch=ops$cnOpts$batch)

+	crlmmOpts <- crlmmOptions(object)$crlmmOpts

+	crlmmOpts$mixtureParams <- mixtureParams

+	crlmmOptions(object)$crlmmOpts <- crlmmOpts

+	pD <- new("AnnotatedDataFrame",

+		  data=sampleStats,

+		  varMetadata=data.frame(labelDescription=colnames(sampleStats)))

+	sampleNames(pD) <- sampleNames(object)

+	phenoData(object) <- pD

+	return(object)

 }

 fitAffySnpMixture56 <- function(S, M, knots, probs=rep(1/3, 3), eps=.01, maxit=10, verbose=FALSE){

@@ -50,13 +50,13 @@ snprma <- function(filenames, mixtureSampleSize=10^5,

 	cnProbes <- get("cnProbes")

 	nr <- length(pnsa) + nrow(cnProbes)

 	nc <- length(filenames)

-	dns <- list(c(gns, rownames(cnProbes)), sns)  

+	dns <- list(c(gns, rownames(cnProbes)), sns)

 	if(isPackageLoaded("ff")){

 		message("ff package is loaded.  Creating ff objects for storing normalized and summarized intensities")

-		A <- ff(dim=c(nr, nc), vmode="integer", finalizer="close", dimnames=dns,

-			overwrite=TRUE)

-		B <- ff(dim=c(nr, nc), vmode="integer", finalizer="close", dimnames=dns,

-			overwrite=TRUE)

+		A <- initializeBigMatrix(nr, nc)

+		B <- initializeBigMatrix(nr, nc)

+		open(A)

+		open(B)

 	}  else {

 		A <- matrix(0L, nr, nc, dimnames=dns)

 		B <- matrix(0L, nr, nc, dimnames=dns)

@@ -98,7 +98,8 @@ snprma <- function(filenames, mixtureSampleSize=10^5,

 	if (!fitMixture) SNR <- mixtureParams <- NA

 	save(A, file=AFile)

 	save(B, file=BFile)

-	return(list(sns=sns, gns=gns, SNR=SNR, SKW=SKW, mixtureParams=mixtureParams, cdfName=cdfName))

+	return(list(sns=sns, gns=gns, SNR=SNR, SKW=SKW, mixtureParams=mixtureParams, cdfName=cdfName,

+		    cnnames=rownames(cnProbes)))

 }

 fitAffySnpMixture56 <- function(S, M, knots, probs=rep(1/3, 3), eps=.01, maxit=10, verbose=FALSE){

@@ -1,85 +1,104 @@

-snprma <- function(filenames, mixtureSampleSize=10^5, fitMixture=FALSE, eps=0.1, verbose=TRUE, seed=1, cdfName, sns){

-  if (missing(sns)) sns <- basename(filenames)

-  ##ADD CHECK TO SEE IF LOADED

-  if (missing(cdfName))

-    cdfName <- read.celfile.header(filenames[1])$cdfName

-  ##  stuffDir <- changeToCrlmmAnnotationName(cdfName)

-  pkgname <- getCrlmmAnnotationName(cdfName)

-  if(!require(pkgname, character.only=TRUE, quietly=!verbose)){

-    suggCall <- paste("library(", pkgname, ", lib.loc='/Altern/Lib/Loc')", sep="")

-    msg <- paste("If", pkgname, "is installed on an alternative location, please load it manually by using", suggCall)

-    message(strwrap(msg))

-    stop("Package ", pkgname, " could not be found.")

-    rm(suggCall, msg)

-  }

-  if(verbose) message("Loading annotations and mixture model parameters.")

-  loader("preprocStuff.rda", .crlmmPkgEnv, pkgname)

-  loader("genotypeStuff.rda", .crlmmPkgEnv, pkgname)

-  loader("mixtureStuff.rda", .crlmmPkgEnv, pkgname)

-  autosomeIndex <- getVarInEnv("autosomeIndex")

-  pnsa <- getVarInEnv("pnsa")

-  pnsb <- getVarInEnv("pnsb")

-  fid <- getVarInEnv("fid")

-  reference <- getVarInEnv("reference")

-  aIndex <- getVarInEnv("aIndex")

-  bIndex <- getVarInEnv("bIndex")

-  SMEDIAN <- getVarInEnv("SMEDIAN")

-  theKnots <- getVarInEnv("theKnots")

-  gns <- getVarInEnv("gns")

+snprma <- function(filenames, mixtureSampleSize=10^5,

+		   fitMixture=FALSE, eps=0.1, verbose=TRUE,

+		   seed=1, cdfName, sns, AFile="A.rda", BFile="B.rda"){

+	if (missing(sns)) sns <- basename(filenames)

+	##ADD CHECK TO SEE IF LOADED

+	if (missing(cdfName))

+		cdfName <- read.celfile.header(filenames[1])$cdfName

+	##  stuffDir <- changeToCrlmmAnnotationName(cdfName)

+	pkgname <- getCrlmmAnnotationName(cdfName)

+	if(!require(pkgname, character.only=TRUE, quietly=!verbose)){

+		suggCall <- paste("library(", pkgname, ", lib.loc='/Altern/Lib/Loc')", sep="")

+		msg <- paste("If", pkgname, "is installed on an alternative location, please load it manually by using", suggCall)

+		message(strwrap(msg))

+		stop("Package ", pkgname, " could not be found.")

+		rm(suggCall, msg)

+	}

+	if(verbose) message("Loading annotations and mixture model parameters.")

+	loader("preprocStuff.rda", .crlmmPkgEnv, pkgname)

+	loader("genotypeStuff.rda", .crlmmPkgEnv, pkgname)

+	loader("mixtureStuff.rda", .crlmmPkgEnv, pkgname)

+	autosomeIndex <- getVarInEnv("autosomeIndex")

+	pnsa <- getVarInEnv("pnsa")

+	pnsb <- getVarInEnv("pnsb")

+	fid <- getVarInEnv("fid")

+	reference <- getVarInEnv("reference")

+	aIndex <- getVarInEnv("aIndex")

+	bIndex <- getVarInEnv("bIndex")

+	SMEDIAN <- getVarInEnv("SMEDIAN")

+	theKnots <- getVarInEnv("theKnots")

+	gns <- getVarInEnv("gns")

-  ##We will read each cel file, summarize, and run EM one by one

-  ##We will save parameters of EM to use later

-  mixtureParams <- matrix(0, 4, length(filenames))

-  SNR <- vector("numeric", length(filenames))

-  SKW <- vector("numeric", length(filenames))

+	##We will read each cel file, summarize, and run EM one by one

+	##We will save parameters of EM to use later

+	mixtureParams <- matrix(0, 4, length(filenames))

+	SNR <- vector("numeric", length(filenames))

+	SKW <- vector("numeric", length(filenames))

-  ## This is the sample for the fitting of splines

-  ## BC: I like better the idea of the user passing the seed,

-  ##     because this might intefere with other analyses

-  ##     (like what happened to GCRMA)

-  set.seed(seed)

-  idx <- sort(sample(autosomeIndex, mixtureSampleSize))

-  ##S will hold (A+B)/2 and M will hold A-B

-  ##NOTE: We actually dont need to save S. Only for pics etc...

-  ##f is the correction. we save to avoid recomputing

-  A <- matrix(as.integer(0), length(pnsa), length(filenames))

-  B <- matrix(as.integer(0), length(pnsb), length(filenames))