@@ -28,3 +28,27 @@ BioC data packages

 * Loaded data in cnrma-functions.R to an environment and extracted from there

   so we can get rid of the NOTES complaining about 'no visible bindings'.

+

+2009-03-18 Rob Scharpf <rscharpf@jhsph.edu> - committed version 1.0.57

+* Added steps object to instantiateObjects to speed up debugging

+* Reformulated the regression for chromosome X

+   -  estimate cross-hybridization using chromosome X

+* Take into account pseudo-autosomal regions on chr X for copy number estimation

+* Replaced get() and assign() with '<-' operations to improve readability

+* Function to compute posterior means of copy number estimates

+

+2009-03-19 Rob Scharpf - committed version 1.0.59

+* Requires genomewidesnp6Crlmm version 1.0.1 or greater

+

+2009-03-25 Rob Scharpf - committed version 1.0.60

+* simplified some of the preliminary steps for the computeCopynumber function

+** crlmm output subset within the body of the function

+* modified vignette -- store results in an oligoSnpSet object (for now)

+* added function to extract the date from the celfile headers (celDates)

+

+

+

+

+

+

+

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

 Package: crlmm

 Type: Package

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

-Version: 1.0.59

+Version: 1.0.60

 Date: 2008-12-28

 Author: Rafael A Irizarry, Benilton S Carvalho <bcarvalh@jhsph.edu>, Robert Scharpf <rscharpf@jhsph.edu>

 Maintainer: Benilton S Carvalho <bcarvalh@jhsph.edu>, Robert Scharpf <rscharpf@jhsph.edu>

@@ -1,5 +1,5 @@

 useDynLib("crlmm", .registration=TRUE)

-export("crlmm", "list.celfiles", "computeCopynumber", "cnrma")

+export("crlmm", "list.celfiles", "computeCopynumber", "cnrma", "celDates")

 exportMethods("list2crlmmSet", "calls", "confs")

 exportClasses("crlmmSet")

 import(Biobase)

@@ -22,21 +22,34 @@ rowCors <- function(x, y, ...){

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

 }

-nuphiAllele <- function(allele, Ystar, W, envir, p){

-	Ns <- get("Ns", envir)

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

-	chrom <- get("chrom", envir)

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

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

+	X <- matrix(x, nrow=nrow)

+	XTX <- crossprod(X)

+	solve(XTX)

+}

+nuphiAllele <- function(p, allele, Ystar, W, envir){

+	Ns <- envir[["Ns"]]

+	CHR <- envir[["chrom"]]

+	nuA <- envir[["nuA"]]

+	nuB <- envir[["nuB"]]

+	nuA.se <- envir[["nuA.se"]]

+	nuB.se <- envir[["nuB.se"]]

+	phiA <- envir[["phiA"]]

+	phiB <- envir[["phiB"]]

+	phiA.se <- envir[["phiA.se"]]

+	phiB.se <- envir[["phiB.se"]]

+	if(CHR==23){

+		phiAx <- envir[["phiAx"]]

+		phiBx <- envir[["phiBx"]]

+	}

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

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

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

-	if(chrom == 23){

-		gender <- get("gender", envir)

-		if(allele == "A" & gender == "female")

-			X <- cbind(1, 2:0)

-		if(allele == "B" & gender == "female")

-			X <- cbind(1, 0:2)

-		if(allele == "A" & gender == "male")

-			X <- cbind(1, 1:0)

-		if(allele == "B" & gender == "male")

-			X <- cbind(1, 0:1)

+	if(CHR == 23){

+		gender <- envir[["gender"]]

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

+		if(allele == "B") X <- cbind(1, c(0, 1, 0, 1, 2), c(1, 0, 2, 1, 0))			

 	} else {##autosome

 		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		

@@ -47,40 +60,57 @@ nuphiAllele <- function(allele, Ystar, W, envir, p){

 		stop("Inf values in W or V")

 	}

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

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

-	##Not instant

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

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

-		X <- matrix(x, nrow=nrow)

-		XTX <- crossprod(X)

-		solve(XTX)

-	}

-	##a little slow

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

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

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

+	if(CHR == 23){

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

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

+	} else{

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

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

+	}

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

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

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

-		ses[, i] <- sqrt(diag(matrix(IXTX[, i], 2, 2) * ssr))

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

+	}

+	if(allele == "A"){

+		nuA[complete, p] <- betahat[1, ]

+		phiA[complete, p] <- betahat[2, ]

+		nuA.se[complete, p] <- ses[1, ]

+		phiA.se[complete, p] <- ses[2, ]

+		envir[["nuA"]] <- nuA

+		envir[["phiA"]] <- phiA

+		envir[["nuA.se"]] <- nuA.se

+		envir[["phiA.se"]] <- phiA.se

+		if(CHR == 23){

+			phiAx[complete, p] <- betahat[3, ]

+			envir[["phiAx"]] <- phiAx

+		}

+	}

+	if(allele=="B"){

+		nuB[complete, p] <- betahat[1, ]

+		phiB[complete, p] <- betahat[2, ]

+		nuB.se[complete, p] <- ses[1, ]

+		phiB.se[complete, p] <- ses[2, ]

+		envir[["nuB"]] <- nuB

+		envir[["phiB"]] <- phiB

+		envir[["nuB.se"]] <- nuB.se

+		envir[["phiB.se"]] <- phiB.se

+		if(CHR == 23){

+			phiBx[complete, p] <- betahat[3, ]

+			envir[["phiBx"]] <- phiBx

+		}

 	}

-	nu <- betahat[1, ]

-	phi <- betahat[2, ]

-	nu.se <- ses[1,]

-	phi.se <- ses[2,]

-	##dimnames(nu) <- dimnames(phi) <- dimnames(nu.ses) <- dimnames(phi.ses)

-	assign("nu", nu, envir=envir)

-	assign("phi", phi, envir=envir)

-	assign("nu.se", nu.se, envir=envir)

-	assign("phi.se", phi.se, envir=envir)	

 }

-celDatesFrom <- function(celfiles, path=""){

+celDates <- function(celfiles){

+	if(!all(file.exists(celfiles))) stop("1 or more cel file does not exist")

 	celdates <- vector("character", length(celfiles))

 	celtimes <- vector("character", length(celfiles))

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

 		if(i %% 100 == 0) cat(".")

-		tmp <- read.celfile.header(file.path(path, celfiles[i]), info="full")$DatHeader

+		tmp <- read.celfile.header(celfiles[i], info="full")$DatHeader

 		tmp <- strsplit(tmp, "\ +")

 		celdates[i] <- tmp[[1]][6]

 		celtimes[i] <- tmp[[1]][7]

@@ -91,14 +121,13 @@ celDatesFrom <- function(celfiles, path=""){

 }

 cnrma <- function(filenames, sns, cdfName, seed=1, verbose=FALSE){

-  ## BC: 03/14/09

-  ## getting pkgname from cdfName, in the future this might be useful

-  ## as the method might be implemented for other platforms

-  pkgname <- getCrlmmAnnotationName(cdfName)

-

+	## BC: 03/14/09

+	## getting pkgname from cdfName, in the future this might be useful

+	## as the method might be implemented for other platforms

+	pkgname <- getCrlmmAnnotationName(cdfName)

 	require(pkgname, character.only=TRUE) || stop("Package ", pkgname, " not available")

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

-  ## Loading data in .crlmmPkgEnv and extracting from there

+	## Loading data in .crlmmPkgEnv and extracting from there

 	data("npProbesFid", package=pkgname, envir=.crlmmPkgEnv)

 	fid <- getVarInEnv("npProbesFid")

 	gc()

@@ -113,7 +142,7 @@ cnrma <- function(filenames, sns, cdfName, seed=1, verbose=FALSE){

 	}

 	##load reference distribution obtained from hapmap

 	data(list="1m_reference_cn", package="genomewidesnp6Crlmm", envir=.crlmmPkgEnv)

-  reference <- getVarInEnv("reference")

+	reference <- getVarInEnv("reference")

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

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

 		x <- log2(y[idx2])

@@ -137,412 +166,471 @@ getFlags <- function(phi.thr, envir){

 	phiB <- get("phiB", envir)

 	negativeNus <- nuA < 1 | nuB < 1

-	negativeNus <- negativeNus > 0

-

 	negativePhis <- phiA < phi.thr | phiB < phi.thr

-	negativePhis <- negativePhis > 0

 	negativeCoef <- negativeNus | negativePhis

 	notfinitePhi <- !is.finite(phiA) | !is.finite(phiB)

-	notfinitePhi <- notfinitePhi > 0

-	flags <- negativeCoef | notfinitePhi  

+	flags <- negativeCoef | notfinitePhi

+	return(flags)

 }

 goodSnps <- function(phi.thr, envir, fewAA=20, fewBB=20){

 	Ns <- get("Ns", envir)

 	flags <- getFlags(phi.thr=phi.thr, envir)

-	fewAA <- Ns[, , 1] < fewAA

-	fewBB <- Ns[, , 3] < fewBB

+	fewAA <- Ns[, , "AA"] < fewAA

+	fewBB <- Ns[, , "BB"] < fewBB

 	flagsA <- flags | fewAA

 	flagsB <- flags | fewBB

 	flags <- list(A=flagsA, B=flagsB)

 	return(flags)

 }

-instantiateObjects <- function(calls, NP, plate, envir, chrom){

-	if(missing(chrom)) stop("must specify chrom")

-	assign("chrom", chrom, envir)

+instantiateObjects <- function(calls, NP, plate, envir, chrom, A, B,

+			       gender, SNRmin=5, SNR){

+	envir[["chrom"]] <- chrom

+	CHR_INDEX <- paste(chrom, "index", sep="")

+	data(list=CHR_INDEX, package="genomewidesnp6Crlmm")	

+	A <- A[index[[1]], SNR > SNRmin]

+	B <- B[index[[1]], SNR > SNRmin]

+	calls <- calls[index[[1]], SNR > SNRmin]

+	conf <- conf[index[[1]], SNR > SNRmin]

+	NP <- NP[index[[2]], SNR > SNRmin]

+	plate <- plate[SNR > SNRmin]

+	uplate <- unique(plate)

+	SNR <- SNR[SNR > SNRmin]

+	

+	envir[["uplate"]] <- uplate

+	envir[["plate"]] <- plate	

+	envir[["NP"]] <- NP

+	envir[["A"]] <- A

+	envir[["B"]] <- B

+	envir[["calls"]] <- calls

+	envir[["conf"]] <- conf

 	snps <- rownames(calls)

 	cnvs <- rownames(NP)

 	sns <- basename(colnames(calls))

 	stopifnot(identical(colnames(calls), colnames(NP)))

-	assign("sns", sns, envir)	

-	assign("snps", snps, envir)

-	assign("cnvs", cnvs, envir)	

+	envir[["sns"]] <- sns

+	envir[["snps"]] <- snps

+	envir[["cnvs"]] <- cnvs

+	if(chrom == 23){

+		if(is.null(gender)){

+			message("Estimating gender")

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

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

+			##gender <- getGender(res)

+			gender[gender==2] <- "female"

+			gender[gender=="1"] <- "male"

+			envir[["gender"]] <- gender

+		} else envir[["gender"]] <- gender

+		phiAx <- matrix(NA, nrow(calls), length(uplate))

+		envir[["phiAx"]] <- phiAx

+		envir[["phiBx"]] <- phiAx

+	}

 	CA <- CB <- matrix(NA, nrow(calls), ncol(calls))

-	assign("plate", plate, envir)

-	uplate <- unique(plate)	

-	Ns <- array(NA, dim=c(nrow(calls), length(uplate), 3))

-	dimnames(Ns)[[3]] <- c("AA", "AB", "BB")

-	muA.AA <- matrix(NA, nrow=nrow(calls), length(uplate))

-	##Sigma.AA <- array(NA, dim=c(nrow(calls), 3, length(uplate)))

-	##dimnames(Sigma.AA)[2:3] <- list(c("varA", "varB", "cov"), uplate)

+	envir[["CA"]] <- CA

+	envir[["CB"]] <- CB

+	

+	Ns <- array(NA, dim=c(nrow(calls), length(uplate), 5))

+	dimnames(Ns)[[3]] <- c("A", "B", "AA", "AB", "BB")

+

+	envir[["Ns"]] <- envir[["muB"]] <- envir[["muA"]] <- Ns

+	envir[["vB"]] <- envir[["vA"]] <- Ns

-	NP.CT <- matrix(NA, nrow(NP), ncol(NP))

-	NP.sds <- matrix(NA, nrow(NP), ncol(NP))

-	assign("NP.CT", NP.CT, envir)

-	assign("NP.sds", NP.sds, envir)

-	nus <- matrix(NA, nrow(NP), length(uplate))

-	assign("nus", nus, envir=envir)  

-	assign("phis", nus, envir=envir)

+	CT.sds <- CT <- matrix(NA, nrow(NP), length(sns))

+	##NP.CT <- matrix(NA, nrow(NP), ncol(NP))

+	##NP.sds <- matrix(NA, nrow(NP), ncol(NP))

+	envir[["CT"]] <- CT

+	envir[["CT.sds"]] <- CT.sds

+

+	##assign("NP.CT", NP.CT, envir)

+	##assign("NP.sds", NP.sds, envir)

+	nuT <- matrix(NA, nrow(NP), length(uplate))

+	phiT <- nuT

+	envir[["nuT"]] <- nuT

+	envir[["phiT"]] <- phiT

+	##assign("nus", nus, envir=envir)  

+	##assign("phis", nus, envir=envir)

 	plates.completed <- rep(FALSE, length(uplate))

-	assign("plates.completed", plates.completed, envir)

-	

-	fit.variance <- NULL

-	assign("fit.variance", fit.variance, envir=envir)

-	npflags <- snpflags <- vector("list", length(uplate))	

-	assign("snpflags", snpflags, envir=envir)

-	assign("npflags", npflags, envir=envir)

-	

-	assign("Ns", Ns, envir=envir)		

-	assign("uplate", uplate, envir=envir)	

-	assign("muA.AA", muA.AA, envir=envir)

-	assign("muA.AB", muA.AA, envir=envir)

-	assign("muA.BB", muA.AA, envir=envir)

-	assign("muB.AA", muA.AA, envir=envir)

-	assign("muB.AB", muA.AA, envir=envir)

-	assign("muB.BB", muA.AA, envir=envir)

-	

-	assign("tau2A", muA.AA, envir=envir)

-	assign("tau2B", muA.AA, envir=envir)

-	assign("sig2A", muA.AA, envir=envir)

-	assign("sig2B", muA.AA, envir=envir)

-	assign("corr", muA.AA, envir=envir)

-	assign("corrA.BB", muA.AA, envir=envir)

-	assign("corrB.AA", muA.AA, envir=envir)			

-	

-	assign("nuA", muA.AA, envir=envir)

-	assign("nuB", muA.AA, envir=envir)

-	assign("phiA", muA.AA, envir=envir)

-	assign("phiB", muA.AA, envir=envir)

-	assign("nuA.se", muA.AA, envir=envir)

-	assign("nuB.se", muA.AA, envir=envir)

-	assign("phiA.se", muA.AA, envir=envir)

-	assign("phiB.se", muA.AA, envir=envir)

-	assign("sd.CT", CA, envir=envir)

-	assign("CA", CA, envir=envir)

-	assign("CB", CB, envir=envir)

+	envir[["plates.completed"]] <- plates.completed

+

+	steps <- rep(FALSE, 4)

+	names(steps) <- c("suffStats", "coef", "snp-cn", "np-cn")

+	envir[["steps"]] <- steps

+

+	snpflags <- matrix(FALSE, length(snps), length(uplate))

+	npflags <- matrix(FALSE, length(cnvs), length(uplate))

+	##assign("snpflags", snpflags, envir=envir)

+	##assign("npflags", npflags, envir=envir)

+	envir[["snpflags"]] <- snpflags

+	envir[["npflags"]] <- npflags

+

+	tau2A <- matrix(NA, nrow(calls), length(uplate))

+	envir[["tau2A"]] <- tau2A

+	envir[["tau2B"]] <- tau2A

+	envir[["sig2A"]] <- tau2A

+	envir[["sig2B"]] <- tau2A

+	sig2T <- matrix(NA, nrow(NP), ncol(NP))

+	envir[["sig2T"]] <- sig2T

+	envir[["corr"]] <- tau2A

+	envir[["corrA.BB"]] <- tau2A

+	envir[["corrB.AA"]] <- tau2A

+	envir[["nuB"]] <- envir[["nuA"]] <- tau2A

+	envir[["phiB"]] <- envir[["phiA"]] <- tau2A

+	envir[["nuB.se"]] <- envir[["nuA.se"]] <- tau2A

+	envir[["phiB.se"]] <- envir[["phiA.se"]] <- tau2A

+	normal <- matrix(TRUE, nrow(A), ncol(A))

+	normalNP <- matrix(TRUE, nrow(NP), ncol(NP))	

+	envir[["normal"]] <- normal

+	envir[["normalNP"]] <- normalNP

 }

-computeCopynumber <- function(A,

+computeCopynumber <- function(chrom,

+			      A,

 			      B,

 			      calls,

 			      conf,

 			      NP,

 			      plate,

-			      fit.variance=NULL,

 			      MIN.OBS=1,

 			      envir,

-			      chrom, P, DF.PRIOR=50, CONF.THR=0.99,

-			      trim=0, upperTail=TRUE,

+			      P,

+			      DF.PRIOR=50,

+			      CONF.THR=0.99,

 			      bias.adj=FALSE,

-			      priorProb, ...){

-	if(length(ls(envir)) == 0) instantiateObjects(calls, NP, plate, envir, chrom)

+			      priorProb,

+			      gender=NULL,

+			      SNR,

+			      SNRmin=5, seed=123, verbose=TRUE, ...){

+	set.seed(seed)

+	if(missing(chrom)) stop("must specify chromosome")

+	if(length(ls(envir)) == 0) {

+		instantiateObjects(calls=calls, NP=NP, plate=plate,

+				   envir=envir, chrom=chrom, A=A, B=B,

+				   gender=gender, SNR=SNR, SNRmin=SNRmin)

+	}

+	plate <- envir[["plate"]]

+	uplate <- envir[["uplate"]]	

+	calls <- envir[["calls"]]

+	A <- envir[["A"]]

+	B <- envir[["B"]]

+	conf <- envir[["conf"]]

+	NP <- envir[["NP"]]

 	if(bias.adj){

 		##assign uniform priors for total copy number states

 		if(missing(priorProb)) priorProb <- rep(1/4, 4)

+		envir[["steps"]] <- rep(FALSE, 4)

 	}

 	##will be updating these objects

-	uplate <- get("uplate", envir)

 	message("Sufficient statistics")

 	if(missing(P)) P <- seq(along=uplate)

-	for(p in P){

-		cat(".")

-		if(sum(plate == uplate[p]) < 10) next()

-		oneBatch(plateIndex=p,

-			 G=calls[, plate==uplate[p]],

-			 A=A[, plate==uplate[p]],

-			 B=B[, plate==uplate[p]],

-			 conf=conf[, plate==uplate[p]],

-			 CONF.THR=CONF.THR,

-			 envir=envir,

-			 MIN.OBS=MIN.OBS,

-			 DF.PRIOR=DF.PRIOR,

-			 trim=trim, upperTail=upperTail,

-			 bias.adj=bias.adj)

+	steps <- envir[["steps"]]

+	if(!steps[1]){

+		for(p in P){

+			cat(".")

+			if(sum(plate == uplate[p]) < 10) next()

+			J <- plate==uplate[p]

+			oneBatch(plateIndex=p,

+				 A=A[, J],

+				 B=B[, J],

+				 calls=calls[, J],

+				 conf=conf[, J],

+				 gender=NULL,

+				 NP[, J],

+				 plate[J],

+				 MIN.OBS=1,

+				 envir=envir,

+				 DF.PRIOR=DF.PRIOR,

+				 CONF.THR=CONF.THR,

+				 bias.adj=bias.adj,

+				 priorProb=priorProb,...)

+		}

+		steps[1] <- TRUE

+		envir[["steps"]] <- steps

 	}

-	message("\nEstimating coefficients")	

-	for(p in P){

-		cat(".")

-		coefs(plateIndex=p, conf=conf[, plate==uplate[p]],

-		      envir=envir, CONF.THR=CONF.THR, MIN.OBS=MIN.OBS)

+	if(!steps[2]){

+		message("\nEstimating coefficients")	

+		for(p in P){

+			cat(".")

+			coefs(plateIndex=p, conf=conf[, plate==uplate[p]],

+			      envir=envir, CONF.THR=CONF.THR, MIN.OBS=MIN.OBS)

+		}

+		steps[2] <- TRUE

+		envir[["steps"]] <- steps		

 	}

-	message("\nAllele specific copy number")	

-	for(p in P){

-		cat(".")

-		polymorphic(plateIndex=p,

-			    A=A[, plate==uplate[p]],

-			    B=B[, plate==uplate[p]],

-			    envir=envir)

+	if(!steps[3]){

+		message("\nAllele specific copy number")	

+		for(p in P){

+			cat(".")

+			polymorphic(plateIndex=p,

+				    A=A[, plate==uplate[p]],

+				    B=B[, plate==uplate[p]],

+				    envir=envir)

+		}

+		steps[3] <- TRUE

+		envir[["snp-cn"]] <- steps						

 	}

-	message("\nCopy number for nonpolymorphic probes...")	

-	for(p in P){

-		cat(".")

-		nonpolymorphic(plateIndex=p,

-			       NP=NP[, plate==uplate[p]],

-			       envir=envir)

+	if(!steps[4]){

+		message("\nCopy number for nonpolymorphic probes...")	

+		for(p in P){

+			cat(".")

+			nonpolymorphic(plateIndex=p,

+				       NP=NP[, plate==uplate[p]],

+				       envir=envir)

+		}

+		steps[4] <- TRUE

+		envir[["np-cn"]] <- steps

 	}

-##	snpflags <- get("snpflags", envir)

-##	npflags <- get("npflags", envir)

-##	flags <- sapply(snpflags, length)

-##	flags.np <- sapply(npflags, length)

-##	if(any(flags > 0) | any(flags.np > 0))

-##		warning("some SNPs were flagged -- possible NAs.  Check the indices in snpflags and npflags")

 }

-nonpolymorphic <- function(plateIndex, NP, envir){

+nonpolymorphic <- function(plateIndex, NP, envir, CONF.THR=0.99, DF.PRIOR=50){

 	p <- plateIndex

-	plate <- get("plate", envir)

-	uplate <- get("uplate", envir)	

-	plates.completed <- get("plates.completed", envir)

+	CHR <- envir[["chrom"]]

+	plate <- envir[["plate"]]

+	uplate <- envir[["uplate"]]

+	plates.completed <- envir[["plates.completed"]]

 	if(!plates.completed[p]) return()

-	##snpflags <- get("snpflags", envir)

-	##if(is.null(snpflags)){

-	snpflags <- goodSnps(phi.thr=10, envir=envir, fewAA=10, fewBB=10)

-	##assign("snpflags", snpflags, envir=envir)

+	snpflags <- goodSnps(phi.thr=2^6, envir=envir, fewAA=10, fewBB=10)

 	flagsA <- snpflags$A[, p]

 	flagsB <- snpflags$B[, p]

 	if(all(flagsA) | all(flagsB)) stop("all snps are flagged")

-	nuA <- get("nuA", envir)

-	nuB <- get("nuB", envir)

-	phiA <- get("phiA", envir)

-	phiB <- get("phiB", envir)

-	uplate <- get("uplate", envir)

-	sns <- get("sns", envir)

-	muA.AA <- get("muA.AA", envir)

-	muA.AB <- get("muA.AB", envir)

-	muA.BB <- get("muA.BB", envir)

-	muB.AA <- get("muB.AA", envir)

-	muB.AB <- get("muB.AB", envir)

-	muB.BB <- get("muB.BB", envir)

-	NP.CT <- get("NP.CT", envir)

-	nus <- get("nus", envir)

-	phis <- get("phis", envir)

-	NP.sds <- get("NP.sds", envir)

-	##fit.variance <- get("fit.variance", envir)

-	NP.CT <- get("NP.CT", envir)

+	nuA <- envir[["nuA"]][, p]

+	nuB <- envir[["nuB"]][, p]

+	phiA <- envir[["phiA"]][, p]

+	phiB <- envir[["phiB"]][, p]

+	uplate <- envir[["uplate"]]

+	sns <- envir[["sns"]]

+	muA <- envir[["muA"]][, p, ]

+	muB <- envir[["muB"]][, p, ]

+	nuT <- envir[["nuT"]]

+	phiT <- envir[["phiT"]]

+	sig2T <- envir[["sig2T"]]

+	A <- envir[["A"]][, plate==uplate[p]]

+	B <- envir[["B"]][, plate==uplate[p]]

+	CA <- envir[["A"]][, plate==uplate[p]]

+	CB <- envir[["B"]][, plate==uplate[p]]

+	if(CHR == 23){

+		phiAx <- envir[["phiAx"]][, p]

+		phiBx <- envir[["phiBx"]][, p]

+	}

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

 	## Train on unflagged SNPs

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

+	##---------------------------------------------------------------------------

+	##Might be best to train using the X chromosome, since for the

+	##X phi and nu have been adjusted for cross-hybridization

 	plateInd <- plate == uplate[p]

-	muA.AAp <- muA.AA[!flagsA, p]

-	muB.BBp <- muB.BB[!flagsB, p]

-		

-	##From SNP data

-	X <- cbind(1, log(c(muA.AAp, muB.BBp)))

-	Y <- log(c(phiA[!flagsA, p], phiB[!flagsB, p]))

-	##Y.nu <- log(c(nuA[!snpflags$A, p], nuB[!snpflags$B, p]))

-	##Y <- cbind(Y.nu, Y.phi)

-	betahat <- solve(crossprod(X), crossprod(X, Y))

-	##Prediction

-	mus <- rowMedians(NP, na.rm=TRUE)

-	X <- cbind(1, log(mus))

-	Yhat <- as.numeric(X %*% betahat)

-	##NP.nus[, p] <- exp(Yhat[, 1])

-	##NP.phis[, p] <- exp(Yhat[, 2])

-	phi <- exp(Yhat)

-	nu <- mus - 2*phi

-	

-	phis[, p] <- as.integer(phi)

-	nus[, p] <- as.integer(nu)

-	##plot(NP.phis[, p], NP.nus[, p], pch=".")

-	CT <- 1/phi*(NP-nu)

-	tmp <- matrix(as.integer(CT*100), nrow(CT), ncol(CT))

-	NP.CT[, plateInd] <- tmp

-	if(FALSE){

-		##should be centered at 2

-		tmp <- rowMeans(NP.CT[, plateInd])

-		hist(tmp/100, breaks=200)

-	}

+	##muA <- muA[!flagsA, p, c("A", "AA")]

+	##muB <- muB[!flagsB, p, c("B", "BB")]

+	muA <- muA[!flagsA, "AA"]

+	muB <- muB[!flagsB, "BB"]

+	X <- cbind(1, log2(c(muA, muB)))

+	Y <- log2(c(phiA[!flagsA], phiB[!flagsB]))

+	if(nrow(X) > 5000) ix <- sample(1:nrow(X), 5000)

+	##X <- X[ix, ]

+	##Y <- Y[ix]

+	betahat <- solve(crossprod(X[ix, ]), crossprod(X[ix, ], Y[ix]))

+##	Yhat <- X%*%betahat

+##	phihat <- 2^Yhat

+##	nuhat <- 2^X[, 2] - 2*phihat

+##	nuAB <- c(nuA[!flagsA], nuB[!flagsB])

+##	plot(log2(nuhat), log2(nuAB), pch=".")

+##	plot(log2(nuhat)-log2(nuAB), pch=".")

+##	hist(log2(nuAB))

+	##plot(Y-Yhat, pch=".")

+	##plot(Y, Yhat, pch=".")

+	##calculate R2

+	if(CHR == 23){

+		cnvs <- envir[["cnvs"]]

+		data(cnProbes, package="genomewidesnp6Crlmm")

+		cnProbes <- cnProbes[match(cnvs, rownames(cnProbes)), ]

+		par <- cnProbes[, "position"] < 2709520 | (cnProbes[, "position"] > 154584237 & cnProbes[, "position"] < 154913754)

+		gender <- envir[["gender"]]

+		mu1 <- rowMedians(NP[, gender=="male"], na.rm=TRUE)

+		mu2 <- rowMedians(NP[, gender=="female"], na.rm=TRUE)

+		mus <- log(cbind(mu1, mu2))

+		X <- cbind(1, mus[, 1])

+		##For build Hg18

+		##http://genome.ucsc.edu/cgi-bin/hgGateway

+		##pseudo-autosomal regions on X

+		##chrX:1-2,709,520 and chrX:154584237-154913754, respectively

+		Yhat1 <- as.numeric(X %*% betahat)

+		X <- cbind(1, mus[, 2])		

+		Yhat2 <- as.numeric(X %*% betahat)

+		phi1 <- exp(Yhat1)

+		phi2 <- exp(Yhat2)

+		nu1 <- exp(mus[, 1]) - phi1

+		nu1[par] <- exp(mus[par, 1]) - 2*phi1[par]

+		nu2 <- exp(mus[, 2]) - 2*phi2

+		CT1 <- 1/phi1*(NP[, gender=="male"]-nu1)

+		CT2 <- 1/phi2*(NP[, gender=="female"]-nu2)

+		CT1 <- matrix(as.integer(100*CT1), nrow(CT1), ncol(CT1))

+		CT2 <- matrix(as.integer(100*CT2), nrow(CT2), ncol(CT2))

+		CT <- envir[["CT"]]

+		CT[, plate==uplate[p] & gender=="male"] <- CT1

+		CT[, plate==uplate[p] & gender=="female"] <- CT2

+		envir[["CT"]] <- CT

+	} else {

+		normalNP <- envir[["normalNP"]]

+		normalNP <- normalNP[, plate==uplate[p]]

+		mus <- rowMedians(NP * normalNP, na.rm=TRUE)

+		crosshyb <- median(muA) - median(mus)

+		X <- cbind(1, log2(mus+crosshyb))

+		##X <- cbind(1, log2(mus))

+		logPhiT <- X %*% betahat

+		phiT[, p] <- 2^(logPhiT)

+		nuT[, p] <- mus - 2*phiT[, p]

+		T <- 1/phiT[, p]*(NP - nuT[, p])

+		CT <- envir[["CT"]]

+		CT[, plate==uplate[p]] <- matrix(as.integer(100*T), nrow(T), ncol(T))

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

-	## For NA SNPs, treat as nonpolymorphic

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

-	CA <- get("CA", envir)

-	CB <- get("CB", envir)	

-	tmpA <- CA[, plate==uplate[p]]

-	tmpB <- CB[, plate==uplate[p]]	

-	indexA <- which(rowSums(is.na(tmpA)) > 0)

-	indexB <- which(rowSums(is.na(tmpB)) > 0)

-	index <- union(indexA, indexB)

-	if(length(index) > 0){

-		npflags <- get("npflags", envir)

-		##warning(paste(length(index), "indices have NAs for the copy number estimates"))		

-		npflags[[p]] <- unique(c(npflags[[p]], index))

-		assign("npflags", npflags, envir)

+		##Variance for prediction region

+		sig2T[, plate==uplate[[p]]] <- rowMAD(log2(NP*normalNP), na.rm=TRUE)^2	

+		envir[["sig2T"]] <- sig2T

+		envir[["CT"]] <- CT

+		envir[["phiT"]] <- phiT

+		envir[["nuT"]] <- nuT

 	}

-	##if(length(index) > 0) browser()

-	

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

-	## this part could stand improvement

-	##  - estimate the uncertainty

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

-	

+	## For NA SNPs, treat as nonpolymorphic

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

-	## Estimate variance for copy number probes

-	## VAR(CT) = var(1/phi * (NP - nu))

-	##         = 1/phi^2 * VAR(NP)

-	##         = 1/phi^2 * f(mu) * sigma

-	## ** Phi is predicted from a regression using the row medians as predictors

-	##   -- this may underestimate the uncertainty

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

-##	log.mus <- log2(mus)

-##	log.var <- log2(rowVars(NP))

-##	f <- predict(fit.variance, newdata=data.frame(log.mus=log.mus))

-##	resid <- log.var - f

-##	sds <- 1/phi*sqrt(2^(predict(fit.variance, newdata=data.frame(log.mus=log.mus+resid))))

-##

-##	robustSD <- function(X) diff(quantile(X, probs=c(0.16, (1-0.16)), na.rm=TRUE))/2

-##	sample.sds <- apply(CT, 2, robustSD)

-##	sample.sds <- sample.sds/median(sample.sds, na.rm=TRUE)

-##	NP.sds[, plate==uplate[p]] <- sds %*% t(sample.sds)

-	assign("phis", phis, envir)

-	assign("nus", nus, envir)

-	assign("NP.CT", NP.CT, envir)

-##	assign("NP.sds", NP.sds, envir)

-##	firstPass.NP <- list(phis=phis, nus=nus, CT=NP.CT, sds=sds,

-##			     gns=gns, sns=sns, bns=bns)

-##	firstPass.NP

 }

-oneBatch <- function(plateIndex, G, A, B, conf, CONF.THR=0.99, MIN.OBS=3, DF.PRIOR, envir, trim, upperTail, bias.adj=FALSE, priorProb, ...){

-	p <- plateIndex

-	plate <- get("plate", envir)

-	Ns <- get("Ns", envir)	

+withinGenotypeMoments <- function(p, A, B, calls, conf, CONF.THR, DF.PRIOR, envir){

+	CHR <- envir[["chrom"]]

+	Ns <- envir[["Ns"]]

+	muA <- envir[["muA"]]

+	muB <- envir[["muB"]]

+	vA <- envir[["vA"]]

+	vB <- envir[["vB"]]

+	plate <- envir[["plate"]]

+	uplate <- envir[["uplate"]]

+	normal <- envir[["normal"]][, plate==uplate[p]]

+	G <- calls; rm(calls); gc()	

+	if(is.null(normal)) normal <- matrix(TRUE, nrow(G), ncol(G))

+

 	highConf <- 1-exp(-conf/1000)

 	highConf <- highConf > CONF.THR

-	AA <- G == 1 & highConf

-	AB <- G == 2 & highConf

-	BB <- G == 3 & highConf

-	Ns[, p, "AA"] <- rowSums(AA)

-	Ns[, p, "AB"] <- rowSums(AB)

-	Ns[, p, "BB"] <- rowSums(BB)

-	assign("Ns", Ns, envir)

-	AA[AA == FALSE] <- NA

-	AB[AB == FALSE] <- NA

-	BB[BB == FALSE] <- NA

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

-	## Sufficient statistics (plate-specific)

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

-	AA.A <- AA*A

-	AB.A <- AB*A

-	BB.A <- BB*A

-	AA.B <- AA*B

-	AB.B <- AB*B

-	BB.B <- BB*B

-	locationAndScale(p=p, AA.A=AA.A, AB.A=AB.A, BB.A=BB.A,

-			 AA.B=AA.B, AB.B=AB.B, BB.B=BB.B,

-			 envir=envir, DF.PRIOR=DF.PRIOR)

-	muA.AA <- get("muA.AA", envir)

-	muA.AB <- get("muA.AB", envir)

-	muA.BB <- get("muA.BB", envir)

-	muB.AA <- get("muB.AA", envir)

-	muB.AB <- get("muB.AB", envir)

-	muB.BB <- get("muB.BB", envir)

-	sigmaA <- get("sigmaA", envir)

-	sigmaB <- get("sigmaB", envir)

-	tau2A <- get("tau2A", envir)

-	tau2B <- get("tau2B", envir)

-	sig2A <- get("sig2A", envir)

-	sig2B <- get("sig2B", envir)

-	corr <- get("corr", envir)

-	corrA.BB <- get("corrA.BB", envir)  ## B allele

-	corrB.AA <- get("corrB.AA", envir)

+	if(CHR == 23){

+		gender <- envir[["gender"]]

+		IX <- matrix(gender, nrow(G), ncol(G), byrow=TRUE)

+		IX <- IX == "female"

+	} else IX <- matrix(TRUE, nrow(G), ncol(G))

+	

+	index <- GT.B <- GT.A <- vector("list", 3)

+	names(index) <- names(GT.B) <- names(GT.A) <- c("AA", "AB", "BB")

+	##--------------------------------------------------

+	##within-genotype sufficient statistics

+	##--------------------------------------------------

+	GT.B <- GT.A <- list()

+	for(j in 1:3){

+		GT <- G==j & highConf & IX & normal

+		Ns[, p, j+2] <- rowSums(GT, na.rm=TRUE)		

+		GT[GT == FALSE] <- NA

+		GT.A[[j]] <- GT*A

+		GT.B[[j]] <- GT*B

+		index[[j]] <- which(Ns[, p, j+2] > 0)

+		

+		muA[, p, j+2] <- rowMedians(GT.A[[j]], na.rm=TRUE)

+		muB[, p, j+2] <- rowMedians(GT.B[[j]], na.rm=TRUE)

+		vA[, p, j+2] <- rowMAD(GT.A[[j]], na.rm=TRUE)

+		vB[, p, j+2] <- rowMAD(GT.B[[j]], na.rm=TRUE)

+

+		##Shrink towards the typical variance

+		DF <- Ns[, p, j+2]-1

+		DF[DF < 1] <- 1

+		v0A <- median(vA[, p, j+2], na.rm=TRUE)

+		v0B <- median(vB[, p, j+2], na.rm=TRUE)

+		vA[, p, j+2] <- (vA[, p, j+2]*DF + v0A*DF.PRIOR)/(DF.PRIOR+DF)

+		vA[is.na(vA[, p, j+2]), p, j+2] <- v0A

+		vB[, p, j+2] <- (vB[, p, j+2]*DF + v0B*DF.PRIOR)/(DF.PRIOR+DF)

+		vB[is.na(vB[, p, j+2]), p, j+2] <- v0B

+	}

+	if(CHR == 23){

+		k <- 1

+		for(j in c(1,3)){

+			GT <- G==j & highConf & !IX 

+			Ns[, p, k] <- rowSums(GT)

+			GT[GT == FALSE] <- NA

+			muA[, p, k] <- rowMedians(GT*A, na.rm=TRUE)

+			muB[, p, k] <- rowMedians(GT*B, na.rm=TRUE)

+			vA[, p, k] <- rowMAD(GT*A, na.rm=TRUE)

+			vB[, p, k] <- rowMAD(GT*B, na.rm=TRUE)

+			

+			DF <- Ns[, p, k]-1

+			DF[DF < 1] <- 1

+			v0A <- median(vA[, p, k], na.rm=TRUE)

+			v0B <- median(vB[, p, k], na.rm=TRUE)

+			vA[, p, k] <- (vA[, p, k]*DF + v0A*DF.PRIOR)/(DF.PRIOR+DF)

+			vA[is.na(vA[, p, k]), p, k] <- v0A

+			vB[, p, k] <- (vB[, p, k]*DF + v0B*DF.PRIOR)/(DF.PRIOR+DF)

+			vB[is.na(vB[, p, k]), p, k] <- v0B			

+			k <- k+1

+		}

+	}

+	envir[["GT.A"]] <- GT.A

+	envir[["GT.B"]] <- GT.B

+	envir[["Ns"]] <- Ns

+	envir[["index"]] <- index

+	envir[["muA"]] <- muA

+	envir[["muB"]] <- muB

+	envir[["vA"]] <- vA

+	envir[["vB"]] <- vB

+}

+	

+

+oneBatch <- function(plateIndex,

+		     A,

+		     B,

+		     calls,

+		     conf,