@@ -516,4 +516,6 @@ function (which expects ff objects and supports parallel processing)

 2010-04-16 B Carvalho committed version 1.5.49

 ** cosmetics - looking for cause of memory spike

+2010-04-24 B Carvalho committed version 1.7.1

+** fixed bug in gender prediction that cause a spike in memory usage

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

 Package: crlmm

 Type: Package

 Title: Genotype Calling (CRLMM) and Copy Number Analysis tool for Affymetrix SNP 5.0 and 6.0 and Illumina arrays.

-Version: 1.7.0

+Version: 1.7.1

 Date: 2010-04-16

 Author: Rafael A Irizarry, Benilton S Carvalho <carvalho@bclab.org>, Robert Scharpf <rscharpf@jhsph.edu>, Matt Ritchie <mritchie@wehi.edu.au>

 Maintainer: Benilton S Carvalho <carvalho@bclab.org>, Robert Scharpf <rscharpf@jhsph.edu>, Matt Ritchie <mritchie@wehi.EDU.AU>

@@ -316,6 +316,27 @@ crlmm2 <- function(filenames, row.names=TRUE, col.names=TRUE,

   return(list2SnpSet(res2, returnParams=returnParams))

 }

+predictGender <- function(cols, theA, theB, XIndex){

+  n <- length(cols)

+  med <- numeric(n)

+  if (n > 0){

+    open(theA)

+    open(theB)

+    for (i in 1:n){

+      vA <- log2(theA[XIndex, cols[i]])

+      vB <- log2(theB[XIndex, cols[i]])

+      med[i] <- median(vA+vB)/2

+      rm(vA, vB)

+    }

+    close(theA)

+    close(theB)

+    rm(theA, theB)

+  }

+  rm(n)

+  gc(verbose=FALSE)

+  med

+}

+

 crlmmGT2 <- function(A, B, SNR, mixtureParams, cdfName, row.names=NULL,

                      col.names=NULL, probs=c(1/3, 1/3, 1/3), DF=6,

                      SNRMin=5, recallMin=10, recallRegMin=1000,

@@ -337,8 +358,8 @@ crlmmGT2 <- function(A, B, SNR, mixtureParams, cdfName, row.names=NULL,

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

   if(verbose) message("Loading annotations.")

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

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

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

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

   ## this is toget rid of the 'no visible binding' notes

   ## variable definitions

   XIndex <- getVarInEnv("XIndex")

@@ -348,12 +369,16 @@ crlmmGT2 <- function(A, B, SNR, mixtureParams, cdfName, row.names=NULL,

   theKnots <- getVarInEnv("theKnots")

   regionInfo <- getVarInEnv("regionInfo")

   params <- getVarInEnv("params")

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

+  rm(obj1, obj2)

   ## IF gender not provide, we predict

   ## FIXME: XIndex may be greater than ocProbesets()

   if(is.null(gender)){

     if(verbose) message("Determining gender.")

-    XMedian <- apply(log2(A[XIndex,, drop=FALSE])+log2(B[XIndex,, drop=FALSE]), 2, median)/2

+##    XMedian <- apply(log2(A[XIndex,, drop=FALSE])+log2(B[XIndex,, drop=FALSE]), 2, median)/2

+    XMedian <- ocLapply(splitIndicesByNode(1:NC), predictGender, theA=A, theB=B, XIndex=XIndex, neededPkgs="crlmm")

+    XMedian <- unlist(XMedian)

     if(sum(SNR[] > SNRMin)==1){

       gender <- which.min(c(abs(XMedian-8.9), abs(XMedian-9.5)))

     }else{

@@ -403,8 +428,8 @@ crlmmGT2 <- function(A, B, SNR, mixtureParams, cdfName, row.names=NULL,

                         sapply(cIndexes, length), SMEDIAN,

                         theKnots, mixtureParams[], DF, probs, 0.025)

-    last <- rSnps[2]

-    rm(snps, rSnps, IndexesBatch, tmpA, tmpB)

+    rm(snps, rSnps, IndexesBatch, tmpA, tmpB, last)

+    gc(verbose=FALSE)

     close(A)

     close(B)

     close(mixtureParams)

@@ -425,7 +450,7 @@ crlmmGT2 <- function(A, B, SNR, mixtureParams, cdfName, row.names=NULL,

   newparams[["centers"]] <- do.call("rbind", lapply(newparamsBatch, "[[", 1))

   newparams[["scales"]] <- do.call("rbind", lapply(newparamsBatch, "[[", 2))

   newparams[["N"]] <- do.call("rbind", lapply(newparamsBatch, "[[", 3))

-  rm(newparamsBatch)

+  rm(newparamsBatch); gc(verbose=FALSE)

   if(verbose) message("Done.")

   if(verbose) message("Estimating recalibration parameters.")

   d <- newparams[["centers"]] - params$centers

@@ -460,6 +485,7 @@ crlmmGT2 <- function(A, B, SNR, mixtureParams, cdfName, row.names=NULL,

       mu <- newparams[["centers"]][i, ]

       j <- which(is.na(mu))

       newparams[["centers"]][i, j] <- c(1, mu[-j], prod(mu[-j]))%*%regParams[, j]

+      rm(mu, j)

     }

     ##remaing NAs are made like originals

@@ -534,8 +560,8 @@ crlmmGT2 <- function(A, B, SNR, mixtureParams, cdfName, row.names=NULL,

                             which(regionInfo[snps, 1]))

     A[snps,] <- tmpA

     B[snps,] <- tmpB

-    last <- rSnps[2]

-    rm(tmpA, tmpB, snps, rSnps, IndexesBatch, ImNull)

+    rm(tmpA, tmpB, snps, rSnps, IndexesBatch, ImNull, last)

+    gc(verbose=FALSE)

     close(A)

     close(B)

     close(mixtureParams)

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

 }