@@ -489,3 +489,10 @@ then readIDAT() should work. Thanks to Pierre Cherel who reported this error.

 ** added parallel/large dataset support to snprma/crlmm

 ** merged changes on .41 with my local changes

+

+2010-04-01 R.Scharpf committed version 1.5.44

+

+** added functions genotype2, cnrma2 for preprocessing and genotyping

+   with crlmm2

+** added crlmmCopynumber2 for parallel support with copy number

+   est. (needs more checking)

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

+Version: 1.5.44

 Date: 2010-02-05

 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>

@@ -12,7 +12,7 @@ importMethodsFrom(Biobase, annotation, "annotation<-",

                   annotatedDataFrameFrom, assayData, "assayData<-",

                   combine, dims, experimentData, "experimentData<-",

                   fData, featureData, "featureData<-", featureNames,

-                  fvarMetadata, fvarLabels, pData, phenoData,

+                  fvarMetadata, fvarLabels, pData, "pData<-", phenoData,

                   "phenoData<-", protocolData, "protocolData<-",

                   pubMedIds, rowMedians, sampleNames, snpCall,

                   snpCallProbability,

@@ -67,7 +67,12 @@ export(crlmm,

        readIdatFiles, 

        snprma,

        snprma2,

-       crlmm2) 

+       crlmm2,

+       genotype2,

+       cnrma2,

+       processCEL2,

+       batch,

+       crlmmCopynumber2)

 export(initializeParamObject)

@@ -1,3 +1,4 @@

+setGeneric("batch", function(object) standardGeneric("batch"))

 setGeneric("getParam", function(object, name, batch) standardGeneric("getParam"))

 setGeneric("cnIndex", function(object) standardGeneric("cnIndex"))

 setGeneric("cnNames", function(object) standardGeneric("cnNames"))

@@ -54,19 +54,20 @@ getFeatureData.Affy <- function(cdfName, copynumber=FALSE){

 	##crlmmOpts$snpRange <- range(snpIndex)

 	##crlmmOpts$npRange <- range(npIndex)

 }

-construct <- function(filenames, cdfName, copynumber=FALSE, sns, verbose=TRUE){

+construct <- function(filenames, cdfName, copynumber=FALSE,

+		      sns, verbose=TRUE){

 	if(verbose) message("Initializing container for assay data elements alleleA, alleleB, call, callProbability")

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

 	protocolData <- getProtocolData.Affy(filenames)

 	featureData <- getFeatureData.Affy(cdfName, copynumber=copynumber)

 	nr <- nrow(featureData); nc <- length(filenames)

 	callSet <- new("SnpSuperSet", 

-			 alleleA=initializeBigMatrix(name="A", nr, nc),

-			 alleleB=initializeBigMatrix(name="B", nr, nc),

-			 call=initializeBigMatrix(name="call", nr, nc),

-			 callProbability=initializeBigMatrix(name="callPr", nr,nc),

-			 featureData=featureData,

-			 annotation=cdfName)

+		       alleleA=initializeBigMatrix(name="A", nr, nc),

+		       alleleB=initializeBigMatrix(name="B", nr, nc),

+		       call=initializeBigMatrix(name="call", nr, nc),

+		       callProbability=initializeBigMatrix(name="callPr", nr,nc),

+		       featureData=featureData,

+		       annotation=cdfName)

 	pd <- data.frame(matrix(NA, nc, 3), row.names=sns)

 	colnames(pd)=c("SKW", "SNR", "gender")

 	phenoData(callSet) <- new("AnnotatedDataFrame", data=pd)

@@ -167,14 +168,557 @@ genotype <- function(filenames,

 			       verbose=verbose,

 			       returnParams=returnParams,

 			       badSNP=badSNP)

-		suppressWarnings(snpCall(callSet)[snp.index, j] <- tmp[["calls"]])

-		suppressWarnings(snpCallProbability(callSet)[snp.index, j] <- tmp[["confs"]])

+		snpCall(callSet)[snp.index, j] <- tmp[["calls"]]

+		snpCallProbability(callSet)[snp.index, j] <- tmp[["confs"]]

 		callSet$gender[j] <- tmp$gender

 		iter <- iter+1

 	}

 	return(callSet)

 }

+genotype2 <- function(filenames,

+		      cdfName,

+		      batch,

+		      mixtureSampleSize=10^5,

+		      eps=0.1,

+		      verbose=TRUE,

+		      seed=1,

+		      sns,

+		      copynumber=FALSE,

+		      probs=rep(1/3, 3),

+		      DF=6,

+		      SNRMin=5,

+		      recallMin=10,

+		      recallRegMin=1000,

+		      gender=NULL,

+		      returnParams=TRUE,

+		      badSNP=0.7){

+	if(!copynumber){

+		callSet <- crlmm2(filenames=filenames,

+				  cdfName=cdfName,

+				  mixtureSampleSize=mixtureSampleSize,

+				  eps=eps,

+				  verbose=verbose,

+				  seed=seed,

+				  sns=sns,

+				  probs=probs,

+				  DF=DF,

+				  SNRMin=SNRMin,

+				  recallMin=recallMin,

+				  recallRegMin=recallRegMin,

+				  gender=gender,

+				  returnParams=returnParams,

+				  badSNP=badSNP)

+		return(callSet)

+	}

+	if(missing(cdfName)) stop("must specify cdfName")

+	if(!isValidCdfName(cdfName)) stop("cdfName not valid.  see validCdfNames")

+	if(missing(batch)){

+		warning("The batch variable is not specified. The scan date of the array will be used as a surrogate for batch.  The batch variable does not affect the preprocessing or genotyping, but is important for copy number estimation.")

+	} else {

+		if(length(batch) != length(filenames))

+			stop("batch variable must be the same length as the filenames")

+	}

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

+	## callSet contains potentially very big matrices

+	callSet <- construct(filenames=filenames,

+			     cdfName=cdfName,

+			     copynumber=TRUE,

+			     sns=sns,

+			     verbose=verbose)

+	if(missing(batch)){

+		protocolData(callSet)$batch <- as.numeric(as.factor(protocolData(callSet)$ScanDate))

+	}

+	if(!missing(batch)) protocolData(callSet)$batch <- batch

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

+	snp.index <- which(isSnp(callSet)==1)

+	snprmaRes <- snprma2(##A=A(callSet),

+			     ##B=B(callSet),

+			       filenames=filenames,

+			       mixtureSampleSize=mixtureSampleSize,

+			       fitMixture=TRUE,

+			       eps=eps,

+			       verbose=verbose,

+			       seed=seed,

+			       cdfName=cdfName,

+			       sns=sns)

+	if(verbose) message("Finished preprocessing.")

+	open(snprmaRes[["A"]])

+	open(snprmaRes[["B"]])

+	open(snprmaRes[["SNR"]])

+	open(snprmaRes[["SKW"]])

+	open(snprmaRes[["mixtureParams"]])

+	if(verbose) message("Updating elements of callSet")

+##	A(callSet) <- snprmaRes[["A"]]

+##	B(callSet) <- snprmaRes[["B"]]

+	bb = ocProbesets()*ncol(A)*8

+	ffrowapply(A(callSet)[i1:i2, ] <- snprmaRes[["A"]][i1:i2, ], X=snprmaRes[["A"]], BATCHBYTES=bb)

+	ffrowapply(B(callSet)[i1:i2, ] <- snprmaRes[["B"]][i1:i2, ], X=snprmaRes[["B"]], BATCHBYTES=bb)

+	##ffrowapply(A(callSet)[i1:i2, ] <- snprmaRes[["A"]][i1:i2, ], X=snprmaRes[["A"]])

+	##ffrowapply(B(callSet)[i1:i2, ] <- snprmaRes[["B"]][i1:i2, ], X=snprmaRes[["B"]])

+	if(verbose) message("Finished updating elements of callSet")

+	stopifnot(identical(featureNames(callSet)[snp.index], snprmaRes$gns))

+	pData(callSet)$SKW <- snprmaRes$SKW

+	pData(callSet)$SNR <- snprmaRes$SNR

+	mixtureParams <- snprmaRes$mixtureParams

+	np.index <- which(isSnp(callSet) == 0)

+	cnrmaRes <- cnrma2(A=A(callSet),

+			   filenames=filenames,

+			   row.names=featureNames(callSet)[np.index],

+			   cdfName=cdfName,

+			   sns=sns,

+			   seed=seed,

+			   verbose=verbose)

+	rm(cnrmaRes); gc()

+	## as.matrix needed when ffdf is used

+	tmp <- crlmmGT2(A=snprmaRes[["A"]],

+			   B=snprmaRes[["B"]],

+			   SNR=snprmaRes[["SNR"]],

+			   mixtureParams=snprmaRes[["mixtureParams"]],

+			   cdfName=cdfName,

+			   row.names=NULL, ##featureNames(callSet),##[snp.index],

+			   col.names=sampleNames(callSet),

+			   probs=probs,

+			   DF=DF,

+			   SNRMin=SNRMin,

+			   recallMin=recallMin,

+			   recallRegMin=recallRegMin,

+			   gender=gender,

+			   verbose=verbose,

+			   returnParams=returnParams,

+			   badSNP=badSNP)

+	open(tmp[["calls"]])

+	open(tmp[["confs"]])

+	##snpCall(callSet) <- tmp[["calls"]]

+	##snpCallProbability(callSet) <- tmp[["confs"]]

+	bb = ocProbesets()*ncol(A)*8

+	ffrowapply(snpCall(callSet)[i1:i2, ] <- tmp[["calls"]][i1:i2, ], X=tmp[["calls"]], BATCHBYTES=bb)

+	ffrowapply(snpCallProbability(callSet)[i1:i2, ] <- tmp[["confs"]][i1:i2, ], X=tmp[["confs"]], BATCHBYTES=bb)

+##	ffrowapply(snpCall(callSet)[i1:i2, ] <- tmp[["calls"]][i1:i2, ], X=tmp[["calls"]])

+##	ffrowapply(snpCallProbability(callSet)[i1:i2, ] <- tmp[["confs"]][i1:i2, ], X=tmp[["confs"]])

+	callSet$gender <- tmp$gender

+	cnSet <- as(callSet, "CNSetLM")

+	return(cnSet)

+}

+

+snprma2RS <- function(A, B,

+		      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_RS, 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_RS <- 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)

+  ##RS

+  iii <- seq(along=pnsa)

+

+  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[iii, k] <- intMedianSummaries(y[aIndex, 1, drop=FALSE], pnsa)

+    B[iii, 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

+}

+

+crlmmGT2_RS<- 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,

+		       gender=NULL, desctrucitve=FALSE, verbose=TRUE,

+		       returnParams=FALSE, badSNP=.7){

+  open(SNR)

+  open(A)

+  open(B)

+  open(mixtureParams)

+  ## expect objects to be ff

+  

+  keepIndex <- which( SNR[] > SNRMin)

+  if(length(keepIndex)==0) stop("No arrays above quality threshold!")

+  

+  NC <- ncol(A)

+  NR <- nrow(A)

+  

+  pkgname <- getCrlmmAnnotationName(cdfName)

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

+

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

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

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

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

+  ## variable definitions

+  XIndex <- getVarInEnv("XIndex")

+  autosomeIndex <- getVarInEnv("autosomeIndex")

+  YIndex <- getVarInEnv("YIndex")

+  SMEDIAN <- getVarInEnv("SMEDIAN")

+  theKnots <- getVarInEnv("theKnots")

+  regionInfo <- getVarInEnv("regionInfo")

+  params <- getVarInEnv("params")

+  ##RS

+  pnsa <- getVarInEnv("pnsa")

+  NR <- length(pnsa)

+  

+  ##IF gender not provide, we predict

+  if(is.null(gender)){

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

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

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

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

+    }else{

+      gender <- kmeans(XMedian, c(min(XMedian[SNR[]>SNRMin]), max(XMedian[SNR[]>SNRMin])))[["cluster"]]

+    }

+  }

+  

+  Indexes <- list(autosomeIndex, XIndex, YIndex)

+  cIndexes <- list(keepIndex, 

+                   keepIndex[which(gender[keepIndex]==2)], 

+                   keepIndex[which(gender[keepIndex]==1)])

+  

+  if(verbose) cat("Calling", NR, "SNPs for recalibration... ")

+

+  ## call C

+  fIndex <- which(gender==2)

+  mIndex <- which(gender==1)

+

+  ## different here

+  ## use gtypeCallerR in batches

+  ##RS

+  snpBatches <- splitIndicesByLength(1:NR, ocProbesets())

+  newparamsBatch <- vector("list", length(snpBatches))

+

+  process1 <- function(idxBatch, snpBatches, autosomeIndex, XIndex,

+                       YIndex, A, B, mixtureParams, fIndex, mIndex,

+                       params, cIndexes, SMEDIAN, theKnots, DF, probs, batchSize){

+    open(A)

+    open(B)

+    open(mixtureParams)

+    snps <- snpBatches[[idxBatch]]

+    rSnps <- range(snps)

+    last <- (idxBatch-1)*batchSize

+    IndexesBatch <- list(autosomeIndex[autosomeIndex %in% snps]-last,

+                         XIndex[XIndex %in% snps]-last,

+                         YIndex[YIndex %in% snps]-last)

+    IndexesBatch <- lapply(IndexesBatch, as.integer)

+    tmpA <- as.matrix(A[snps,])

+    tmpB <- as.matrix(B[snps,])

+    ## newparamsBatch[[idxBatch]]

+

+    tmp <- gtypeCallerR(tmpA, tmpB, fIndex, mIndex,

+                        params[["centers"]][snps,],

+                        params[["scales"]][snps,],

+                        params[["N"]][snps,],

+                        IndexesBatch, cIndexes,

+                        sapply(IndexesBatch, length),

+                        sapply(cIndexes, length), SMEDIAN,

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

+    

+    last <- rSnps[2]

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

+    close(A)

+    close(B)

+    close(mixtureParams)

+    tmp

+  }

+

+  newparamsBatch <- ocLapply(seq(along=snpBatches), process1,

+                             snpBatches=snpBatches,

+                             autosomeIndex=autosomeIndex, XIndex=XIndex,

+                             YIndex=YIndex, A=A, B=B,

+                             mixtureParams=mixtureParams, fIndex=fIndex,

+                             mIndex=mIndex, params=params,

+                             cIndexes=cIndexes, SMEDIAN=SMEDIAN,

+                             theKnots=theKnots, DF=DF, probs=probs,

+                             batchSize=ocProbesets())

+##   last <- 0

+##   for (idxBatch in seq(along=snpBatches)){

+##     snps <- snpBatches[[idxBatch]]

+##     rSnps <- range(snps)

+##     IndexesBatch <- list(autosomeIndex[autosomeIndex %in% snps]-last,

+##                          XIndex[XIndex %in% snps]-last,

+##                          YIndex[YIndex %in% snps]-last)

+##     IndexesBatch <- lapply(IndexesBatch, as.integer)

+##     tmpA <- A[snps,]

+##     tmpB <- B[snps,]

+##     newparamsBatch[[idxBatch]] <- gtypeCallerR(tmpA, tmpB, fIndex, mIndex,

+##                                                params[["centers"]][snps,],

+##                                                params[["scales"]][snps,],

+##                                                params[["N"]][snps,],

+##                                                IndexesBatch, cIndexes,

+##                                                sapply(IndexesBatch, length),

+##                                                sapply(cIndexes, length),

+##                                                SMEDIAN, theKnots,

+##                                                mixtureParams[], DF, probs, 0.025)

+##     last <- rSnps[2]

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

+##   }

+##   rm(last)

+  

+  newparams <- vector("list", 3)

+  names(newparams) <- c("centers", "scales", "N")

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

+  if(verbose) message("Done.")

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

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

+

+  ##regression 

+  Index <- intersect(which(pmin(newparams[["N"]][, 1],

+                                newparams[["N"]][, 2],

+                                newparams[["N"]][, 3]) > recallMin &

+                                !apply(regionInfo, 1, any)),

+                                autosomeIndex)

+  if(length(Index) < recallRegMin){

+    warning("Recalibration not possible. Possible cause: small sample size.")

+    newparams <- params

+    dev <- vector("numeric", nrow(newparams[["centers"]]))

+    SS <- matrix(Inf, 3, 3)

+    DD <- 0

+  }else{

+    data4reg <- as.data.frame(newparams[["centers"]][Index,])

+    names(data4reg) <- c("AA", "AB", "BB")

+    regParams <- cbind(  coef(lm(AA~AB*BB, data=data4reg)),

+                       c(coef(lm(AB~AA+BB, data=data4reg)), 0), 

+                       coef(lm(BB~AA*AB, data=data4reg)))

+    rownames(regParams) <- c("intercept", "X", "Y", "XY")

+    rm(data4reg)

+  

+    minN <- 3

+    newparams[["centers"]][newparams[["N"]] < minN] <- NA

+    Index <- setdiff(which(rowSums(is.na(newparams[["centers"]]))==1), YIndex)

+    if(verbose) cat("Filling out empty centers")

+    for(i in Index){

+      if(verbose) if(i%%10000==0)cat(".")

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

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

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

+    }

+    

+    ##remaing NAs are made like originals

+    if(length(YIndex)>0){

+      noMoveIndex <- union(setdiff(which(rowSums(is.na(newparams[["centers"]]))>0), YIndex), 

+                           YIndex[rowSums(is.na(newparams[["centers"]][YIndex, ])>1)])

+    }

+    snps2ignore <- which(rowSums(is.na(newparams[["centers"]])) > 0)

+    snps2keep <- setdiff(autosomeIndex, snps2ignore)

+    rm(snps2ignore)

+    newparams[["centers"]][is.na(newparams[["centers"]])] <- params[["centers"]][is.na(newparams[["centers"]])]

+    if(verbose) cat("\n")

+  

+    if(verbose) message("Calculating and standardizing size of shift.")

+    GG <- DD <- newparams[["centers"]] - params[["centers"]]

+    DD <- sweep(DD, 2, colMeans(DD[autosomeIndex, ]))

+    SS <- cov(DD[autosomeIndex, ])

+    SSI <- solve(SS)

+    dev <- vector("numeric", nrow(DD))

+    if(length(YIndex)){

+      dev[-YIndex] <- apply(DD[-YIndex, ], 1, function(x) x%*%SSI%*%x)

+      dev[-YIndex] <- 1/sqrt( (2*pi)^3*det(SS))*exp(-0.5*dev[-YIndex])

+      ##Now Y (only two params)

+      SSY <- SS[c(1, 3), c(1, 3)]

+      SSI <- solve(SSY) 

+      dev[YIndex] <- apply(DD[YIndex, c(1, 3)], 1, function(x) x%*%SSI%*%x)

+      dev[YIndex] <- 1/sqrt( (2*pi)^2*det(SSY))*exp(-0.5*dev[YIndex])

+    } else {

+      dev=apply(DD,1,function(x) x%*%SSI%*%x)

+      dev=1/sqrt( (2*pi)^3*det(SS))*exp(-0.5*dev)

+    }

+  }

+    

+  ## BC: must keep SD

+  params[-2] <- newparams[-2]

+  

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

+  if(verbose) cat("Calling", NR, "SNPs... ")

+  ## ###################

+  ## ## MOVE TO C#######

+

+  ## running in batches

+  ## snpBatches <- splitIndicesByLength(1:nrow(A), ocProbesets())

+

+  process2 <- function(idxBatch, snpBatches, autosomeIndex, XIndex,

+                       YIndex, A, B, mixtureParams, fIndex, mIndex,

+                       params, cIndexes, SMEDIAN, theKnots, DF, probs,

+                       regionInfo, batchSize){

+    open(A)

+    open(B)

+    open(mixtureParams)

+    snps <- snpBatches[[idxBatch]]

+    tmpA <- as.matrix(A[snps,])

+    tmpB <- as.matrix(B[snps,])

+    rSnps <- range(snps)

+    last <- (idxBatch-1)*batchSize

+    IndexesBatch <- list(autosomeIndex[autosomeIndex %in% snps]-last,

+                         XIndex[XIndex %in% snps]-last,

+                         YIndex[YIndex %in% snps]-last)

+    IndexesBatch <- lapply(IndexesBatch, as.integer)

+    ImNull <- gtypeCallerR2(tmpA, tmpB, fIndex, mIndex,

+                            params[["centers"]][snps,],

+                            params[["scales"]][snps,],

+                            params[["N"]][snps,],

+                            IndexesBatch, cIndexes,

+                            sapply(IndexesBatch, length),

+                            sapply(cIndexes, length),

+                            SMEDIAN, theKnots, mixtureParams[],

+                            DF, probs, 0.025,

+                            which(regionInfo[snps, 2]),

+                            which(regionInfo[snps, 1]))

+    A[snps,] <- tmpA

+    B[snps,] <- tmpB

+    last <- rSnps[2]

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

+    close(A)

+    close(B)

+    close(mixtureParams)

+  }

+

+  ocLapply(seq(along=snpBatches), process2, snpBatches=snpBatches,

+           autosomeIndex=autosomeIndex, XIndex=XIndex, YIndex=YIndex,

+           A=A, B=B, mixtureParams=mixtureParams, fIndex=fIndex,

+           mIndex=mIndex, params=params, cIndexes=cIndexes,

+           SMEDIAN=SMEDIAN, theKnots=theKnots, DF=DF, probs=probs,

+           regionInfo=regionInfo, batchSize=ocProbesets())

+  

+##   last <- 0

+##   for (idxBatch in seq(along=snpBatches)){

+##     snps <- snpBatches[[idxBatch]]

+##     tmpA <- A[snps,]

+##     tmpB <- B[snps,]

+##     rSnps <- range(snps)

+##     IndexesBatch <- list(autosomeIndex[autosomeIndex %in% snps]-last,

+##                          XIndex[XIndex %in% snps]-last,

+##                          YIndex[YIndex %in% snps]-last)

+##     IndexesBatch <- lapply(IndexesBatch, as.integer)

+##     ImNull <- gtypeCallerR2(tmpA, tmpB, fIndex, mIndex,

+##                             params[["centers"]][snps,],

+##                             params[["scales"]][snps,],

+##                             params[["N"]][snps,],

+##                             IndexesBatch, cIndexes,

+##                             sapply(IndexesBatch, length),

+##                             sapply(cIndexes, length),

+##                             SMEDIAN, theKnots, mixtureParams[],

+##                             DF, probs, 0.025,

+##                             which(regionInfo[snps, 2]),

+##                             which(regionInfo[snps, 1]))

+##     A[snps,] <- tmpA

+##     B[snps,] <- tmpB

+##     last <- rSnps[2]

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

+##   }

+##   

+##   gc(verbose=FALSE)

+  ##  END MOVE TO C#######

+  ## ##################

+  

+  dev <- dev/(dev+1/383)

+  if(!is.null(row.names)){ rownames(A) <- rownames(B) <- names(dev) <- row.names}

+  if(!is.null(col.names)){ colnames(A) <- colnames(B) <- col.names}

+

+  if(length(Index) >= recallRegMin){

+   tmp4batchQC <- DD[autosomeIndex,]/(params[["N"]][autosomeIndex,]+1)

+   tmpSnpQc <- dev[autosomeIndex]

+   SS <- cov(tmp4batchQC[tmpSnpQc < badSNP,])

+   batchQC <- mean(diag(SS))

+  }else{

+    SS <- matrix(0, 3, 3)

+    batchQC <- Inf

+  }

+  

+  if(verbose) message("Done.")

+  if (returnParams){

+    return(list(calls=A, confs=B, SNPQC=dev, batchQC=batchQC, params=params, DD=DD, covDD=SS, gender=gender, pkgname=pkgname))

+  }else{

+    return(list(calls=A, confs=B, SNPQC=dev, batchQC=batchQC, DD=DD, covDD=SS, gender=gender, pkgname=pkgname))

+  }

+}

+

+

 ##---------------------------------------------------------------------------

 ##---------------------------------------------------------------------------

 ## For Illumina

@@ -252,6 +796,7 @@ constructRG <- function(filenames, cdfName, sns, verbose, fileExt, sep, sampleSh

 }

 crlmmIlluminaRS <- function(sampleSheet=NULL,

 			    arrayNames=NULL,

+			    batch,

 			    ids=NULL,

 			    path=".",

 			    arrayInfoColNames=list(barcode="SentrixBarcode_A", position="SentrixPosition_A"),

@@ -266,10 +811,18 @@ crlmmIlluminaRS <- function(sampleSheet=NULL,

 			    seed=1, save.ab=FALSE, snpFile, cnFile,

 			    mixtureSampleSize=10^5, eps=0.1, verbose=TRUE,

 			    cdfName, sns, recallMin=10, recallRegMin=1000,

-			    returnParams=FALSE, badSNP=.7) {

+			    returnParams=FALSE, badSNP=.7,

+			    copynumber=FALSE,

+			    load.it=TRUE) {

 	if(missing(cdfName)) stop("must specify cdfName")

 	if(!isValidCdfName(cdfName)) stop("cdfName not valid.  see validCdfNames")

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

+	if(missing(batch)){

+		warning("The batch variable is not specified. The scan date of the array will be used as a surrogate for batch.  The batch variable does not affect the preprocessing or genotyping, but is important for copy number estimation.")

+	} else {

+		if(length(batch) != length(sns))

+			stop("batch variable must be the same length as the filenames")

+	}	

 	batches <- splitIndicesByLength(seq(along=arrayNames), ocSamples())

 	k <- 1

 	for(j in batches){

@@ -300,13 +853,17 @@ crlmmIlluminaRS <- function(sampleSheet=NULL,

 		##  Here, I'm just using the # of rows returned from the above function

 		if(k == 1){

 			if(verbose) message("Initializing container for alleleA, alleleB, call, callProbability")

-			callSet <- new("SnpSuperSet",

-				       alleleA=initializeBigMatrix(name="A", nr=nrow(res[[1]]), nc=length(sns)),

-				       alleleB=initializeBigMatrix(name="B", nr=nrow(res[[1]]), nc=length(sns)),

-				       call=initializeBigMatrix(name="call", nr=nrow(res[[1]]), nc=length(sns)),

-				       callProbability=initializeBigMatrix(name="callPr", nr=nrow(res[[1]]), nc=length(sns)),

-				       annotation=cdfName)

-			sampleNames(callSet) <- sns

+			load.obj <- loadObject("callSet", load.it)

+			if(!load.obj){

+				callSet <- new("SnpSuperSet",

+					       alleleA=initializeBigMatrix(name="A", nr=nrow(res[[1]]), nc=length(sns)),

+					       alleleB=initializeBigMatrix(name="B", nr=nrow(res[[1]]), nc=length(sns)),

+					       call=initializeBigMatrix(name="call", nr=nrow(res[[1]]), nc=length(sns)),

+					       callProbability=initializeBigMatrix(name="callPr", nr=nrow(res[[1]]), nc=length(sns)),

+					       annotation=cdfName)

+				sampleNames(callSet) <- sns

+				save(callSet, file=file.path(ldPath(), "callSet.rda"))

+			} else load(file.path(ldPath(), "callSet.rda"))

 			phenoData(callSet) <- getPhenoData(sampleSheet=sampleSheet,

 							   arrayNames=sns,

 							   arrayInfoColNames=arrayInfoColNames)

@@ -315,9 +872,21 @@ crlmmIlluminaRS <- function(sampleSheet=NULL,

 			protocolData(callSet) <- new("AnnotatedDataFrame", data=pD)

 			pData(protocolData(callSet))[j, ] <- pData(protocolData)

 			featureNames(callSet) <- res[["gns"]]

-			pData(callSet)$SKW <- rep(NA, length(sns))

-			pData(callSet)$SNR <- rep(NA, length(sns))

+			pData(callSet)$SNR <- initializeBigVector("crlmmSNR-", length(sns), "double")

+			pData(callSet)$SKW <- initializeBigVector("crlmmSKW-", length(sns), "double")

+			##pData(callSet)$SKW <- rep(NA, length(sns))

+			##pData(callSet)$SNR <- rep(NA, length(sns))

 			pData(callSet)$gender <- rep(NA, length(sns))

+			mixtureParams <- initializeBigMatrix("crlmmMixt-", nr=4, nc=ncol(callSet), vmode="double")

+			if(missing(batch)){

+				protocolData(callSet)$batch <- rep(NA, length(sns))

+			} else{

+				protocolData(callSet)$batch <- batch

+			}

+			featureData(callSet) <- addFeatureAnnotation(callSet)

+			open(mixtureParams)

+			open(callSet$SNR)

+			open(callSet$SKW)

 		}

 		if(k > 1 & nrow(res[[1]]) != nrow(callSet)){

 			##RS: I don't understand why the IDATS for the

@@ -325,40 +894,61 @@ crlmmIlluminaRS <- function(sampleSheet=NULL,

 			res[["A"]] <- res[["A"]][res$gns %in% featureNames(callSet), ]

 			res[["B"]] <- res[["B"]][res$gns %in% featureNames(callSet), ]

 		}

+		if(missing(batch)){

+			protocolData(callSet)$batch[j] <- as.numeric(as.factor(protocolData$ScanDate))

+		}

 		## MR: we need to define a snp.index vs np.index

 		snp.index <- match(res$gns, featureNames(callSet))		

-		suppressWarnings(A(callSet)[snp.index, j] <- res[["A"]])

-		suppressWarnings(B(callSet)[snp.index, j] <- res[["B"]])

+		A(callSet)[snp.index, j] <- res[["A"]]

+		B(callSet)[snp.index, j] <- res[["B"]]

 		pData(callSet)$SKW[j] <- res$SKW

 		pData(callSet)$SNR[j] <- res$SNR

-		mixtureParams <- res$mixtureParams

+		mixtureParams[, j] <- res$mixtureParams

 		rm(res); gc()

-		##MR:  edit snp.index

-		##snp.index <- 1:nrow(callSet)

-		tmp <- crlmmGT(A=as.matrix(A(callSet)[, j]),

-			       B=as.matrix(B(callSet)[, j]),

-			       SNR=callSet$SNR[j],

-			       mixtureParams=mixtureParams,

-			       cdfName=annotation(callSet),

-			       row.names=featureNames(callSet),

-			       col.names=sampleNames(callSet)[j],

-			       probs=probs,

-			       DF=DF,

-			       SNRMin=SNRMin,

-			       recallMin=recallMin,

-			       recallRegMin=recallRegMin,

-			       gender=gender,

-			       verbose=verbose,

-			       returnParams=returnParams,

-			       badSNP=badSNP)

-		suppressWarnings(snpCall(callSet)[, j] <- tmp[["calls"]])

-		## MR: many zeros in the conf. scores (?)

-		suppressWarnings(snpCallProbability(callSet)[, j] <- tmp[["confs"]])

-		callSet$gender[j] <- tmp$gender

-		rm(tmp); gc()

 		k <- k+1

 	}

-	return(callSet)

+	tmp <- crlmmGT2(A=A(callSet),

+			B=B(callSet),

+			SNR=callSet$SNR,

+			mixtureParams=mixtureParams,

+			cdfName=annotation(callSet),

+			row.names=featureNames(callSet),

+			col.names=sampleNames(callSet),

+			probs=probs,

+			DF=DF,

+			SNRMin=SNRMin,

+			recallMin=recallMin,

+			recallRegMin=recallRegMin,

+			gender=gender,

+			verbose=verbose,

+			returnParams=returnParams,

+			badSNP=badSNP)

+	open(tmp[["calls"]])

+	open(tmp[["confs"]])

+	snpCall(callSet) <- tmp[["calls"]]

+	## MR: many zeros in the conf. scores (?)

+	snpCallProbability(callSet) <- tmp[["confs"]]

+	callSet$gender <- tmp$gender

+	if(copynumber){

+		load.obj <- loadObject("cnSet", load.it)

+		if(!load.obj){

+			cnSet <- as(callSet, "CNSetLM")

+		} else {

+			load(file.path(ldPath(), "cnSet.rda"))

+			A(cnSet) <- A(callSet)

+			B(cnSet) <- B(callSet)

+			snpCall(cnSet) <- snpCall(callSet)

+			snpCallProbability(cnSet) <- snpCallProbability(callSet)

+			annotation(cnSet) <- annotation(callSet)

+			featureData(cnSet) <- featureData(callSet)

+			protocolData(cnSet) <- protocolData(callSet)

+			phenoData(cnSet) <- phenoData(callSet)

+			experimentData(cnSet) <- experimentData(callSet)

+		}

+	}

+	close(mixtureParams)

+	rm(tmp); gc()

+	return(cnSet)

 }

 ##---------------------------------------------------------------------------

 ##---------------------------------------------------------------------------

@@ -377,6 +967,29 @@ rowMAD <- function(x, y, ...){

 	return(mad)

 }

+shrink <- function(x, Ns, DF.PRIOR){

+	DF <- Ns-1

+	DF[DF < 1] <- 1

+	x.0 <- apply(x, 2, median, na.rm=TRUE)

+	x <- (x*DF + x.0*DF.PRIOR)/(DF.PRIOR + DF)

+	for(j in 1:ncol(x)) x[is.na(x[, j]), j] <- x.0[j]

+	return(x)

+}

+

+applyByGenotype <- function(x, FUN, G){

+	FUN <- match.fun(FUN)

+	tmp <- matrix(NA, nrow(x), 3)

+	for(j in 1:3){

+		GT <- G==j

+		GT[GT == FALSE] <- NA

+		gt.x <- GT*x

+		tmp[, j] <- FUN(gt.x, na.rm=TRUE)

+	}

+	tmp

+}

+

+

+

 rowCors <- function(x, y, ...){

 	N <- rowSums(!is.na(x))

 	x <- suppressWarnings(log2(x))

@@ -387,6 +1000,15 @@ rowCors <- function(x, y, ...){

 	return(covar/(sd.x*sd.y))

 }

+corByGenotype <- function(A, B, G, Ns, which.cluster=c(1,2,3)[1], DF.PRIOR){

+	x <- A * (G == which.cluster)

+	x[x==0] <- NA

+	y <- B * (G == which.cluster)

+	res <- as.matrix(rowCors(x, y, na.rm=TRUE))

+	cors <- shrink(res, Ns[, which.cluster], DF.PRIOR)

+	cors

+}

+

 generateX <- function(w, X) as.numeric(diag(w) %*% X)

 generateIXTX <- function(x, nrow=3) {

 	X <- matrix(x, nrow=nrow)

@@ -394,6 +1016,61 @@ generateIXTX <- function(x, nrow=3) {

 	solve(XTX)

 }

+nuphiAllele2 <- function(allele, Ystar, W, Ns){

+	complete <- rowSums(is.na(W)) == 0 

+	if(any(!is.finite(W))){## | any(!is.finite(V))){

+		i <- which(rowSums(!is.finite(W)) > 0)

+		stop("Possible zeros in the within-genotype estimates of the spread (vA, vB). ")

+	}

+	NOHET <- mean(Ns[, 2], na.rm=TRUE) < 0.05

+	if(missing(allele)) stop("must specify allele")

+	if(allele == "A") X <- cbind(1, 2:0) else X <- cbind(1, 0:2)

+	if(NOHET) X <- X[-2, ] ##more than 1 X chromosome, but all homozygous		

+	##How to quickly generate Xstar, Xstar = diag(W) %*% X

+	Xstar <- apply(W, 1, generateX, X)

+	IXTX <- apply(Xstar, 2, generateIXTX, nrow=nrow(X))

+	betahat <- matrix(NA, 2, nrow(Ystar))

+	ses <- matrix(NA, 2, nrow(Ystar))

+	for(i in 1:nrow(Ystar)){

+		betahat[, i] <- crossprod(matrix(IXTX[, i], ncol(X), ncol(X)), crossprod(matrix(Xstar[, i], nrow=nrow(X)), Ystar[i, ]))

+		ssr <- sum((Ystar[i, ] - matrix(Xstar[, i], nrow(X), ncol(X)) %*% matrix(betahat[, i], ncol(X), 1))^2)

+		ses[, i] <- sqrt(diag(matrix(IXTX[, i], ncol(X), ncol(X)) * ssr))

+	}

+	nu <- betahat[1, ]

+	phi <- betahat[2, ]

+	return(list(nu, phi))

+}

+

+nuphiAlleleX <- function(allele, Ystar, W, Ns, chrX=FALSE){

+	complete <- rowSums(is.na(W)) == 0 

+	if(any(!is.finite(W))){## | any(!is.finite(V))){