@@ -1,12 +1,12 @@

 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.15.0

+Version: 1.15.12

 Author: Benilton S Carvalho, Robert Scharpf, Matt Ritchie, Ingo Ruczinski, Rafael A Irizarry

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

 Description: Faster implementation of CRLMM specific to SNP 5.0 and 6.0 arrays, as well as a copy number tool specific to 5.0, 6.0, and Illumina platforms

 License: Artistic-2.0

-Depends: R (>= 2.14.0), oligoClasses (>= 1.17.36)

+Depends: R (>= 2.14.0), oligoClasses (>= 1.19.17)

 Imports: methods,

          Biobase (>= 2.15.4),

          BiocGenerics,

@@ -20,7 +20,8 @@ Imports: methods,

          SNPchip,

          utils,

 	 lattice,

-	 ff

+	 ff,

+	 foreach

 Suggests: hapmapsnp6,

           genomewidesnp6Crlmm (>= 1.0.4),

           GGdata,

@@ -53,7 +53,7 @@ importFrom(oligoClasses, celfileDate, chromosome2integer, i2p,

 	   isPackageLoaded,

            ldPath, ocLapply, ocProbesets, ocSamples,

 	   parStatus,

-           splitIndicesByLength, splitIndicesByNode)

+           splitIndicesByLength, splitIndicesByNode, AssayDataList)

 importFrom(preprocessCore, normalize.quantiles,

            normalize.quantiles.use.target)

@@ -64,11 +64,14 @@ importFrom(stats, coef, cov, dnorm, kmeans, lm, mad, median, quantile,

 importFrom(utils, packageDescription, setTxtProgressBar,

            txtProgressBar)

+## foreach

+import(foreach)

+

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

 ## export

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

 exportClasses(PredictionRegion)

-exportMethods(CA, CB, coerce,

+exportMethods(CA, CB,

 	      A, B, corr, nuA, nuB, phiA, phiB,

 	      predictionRegion, posteriorProbability,

 	      tau2, Ns, medians, mads,

@@ -87,4 +90,4 @@ export(crlmm,

        genotype.Illumina,

        crlmmCopynumber2, crlmmCopynumberLD, crlmmCopynumber)

 export(genotypes, totalCopynumber, rawCopynumber, xyplot)

-export(ABpanel, validCEL, celDates)

+export(ABpanel, validCEL, celDates, validCdfNames)

@@ -104,4 +104,4 @@ setGeneric("predictionRegion", function(object, copyNumber=0:4)

 	   standardGeneric("predictionRegion"))

 setGeneric("xyplot", useAsDefault=function(x, data, ...) lattice::xyplot(x, data,...))

 setGeneric("xyplotcrlmm", function(x, data, predictRegion,...) standardGeneric("xyplotcrlmm"))

-setGeneric("calculateRBaf", function(object, batch.name) standardGeneric("calculateRBaf"))

+setGeneric("calculateRBaf", function(object, batch.name, chrom) standardGeneric("calculateRBaf"))

@@ -21,7 +21,7 @@ getProtocolData.Affy <- function(filenames){

 ##		  return(gns)

 ##	  })

-getFeatureData <- function(cdfName, copynumber=FALSE){

+getFeatureData <- function(cdfName, copynumber=FALSE, genome=genome){

 	pkgname <- getCrlmmAnnotationName(cdfName)

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

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

@@ -39,10 +39,15 @@ getFeatureData <- function(cdfName, copynumber=FALSE){

 		pkgname <- paste(pkgname, "Crlmm", sep="")

 	}

 	path <- system.file("extdata", package=pkgname)

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

+	multiple.builds <- length(grep("hg19", list.files(path)) > 0)

+	if(!multiple.builds)

+		load(file.path(path, "snpProbes.rda"))

+	else load(file.path(path, paste("snpProbes_", genome, ".rda", sep="")))

 	snpProbes <- get("snpProbes")

 	if(copynumber){

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

+		if(!multiple.builds)

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

+		else load(file.path(path, paste("cnProbes_", genome, ".rda", sep="")))

 		cnProbes <- get("cnProbes")

 		snpIndex <- seq(along=gns)

 		npIndex <- seq(along=rownames(cnProbes)) + max(snpIndex)

@@ -121,8 +126,8 @@ construct <- function(filenames,

 	return(cnSet)

 }

-constructAffy <- function(filenames, sns, cdfName, batch, verbose=TRUE){

-	stopifnot(!missing(filenames))

+constructAffy <- function(filenames, sns, cdfName, batch, verbose=TRUE, genome){

+	##stopifnot(!missing(filenames))

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

 	stopifnot(!missing(batch))

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

@@ -131,7 +136,6 @@ constructAffy <- function(filenames, sns, cdfName, batch, verbose=TRUE){

 	## appropriate dimension for snps+nps

 	##

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

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

 	if (missing(cdfName))

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

 	pkgname <- getCrlmmAnnotationName(cdfName)

@@ -148,7 +152,7 @@ constructAffy <- function(filenames, sns, cdfName, batch, verbose=TRUE){

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

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

 	genderff <- initializeBigVector("gender", length(filenames), "integer")

-	featureData <- getFeatureData(cdfName, copynumber=TRUE)

+	featureData <- getFeatureData(cdfName, copynumber=TRUE, genome=genome)

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

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

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

@@ -164,18 +168,12 @@ constructAffy <- function(filenames, sns, cdfName, batch, verbose=TRUE){

 		     callProbability=callPr,

 		     featureData=featureData,

 		     annotation=cdfName,

-		     batch=as.character(batch))

-	if(!missing(sns)){

-		sampleNames(cnSet) <- sns

-	} else {

-		sampleNames(cnSet) <- basename(filenames)

-	}

+		     batch=as.character(batch),

+		     genome=genome)

+	sampleNames(cnSet) <- sns

 	protocolData <- getProtocolData.Affy(filenames)

 	rownames(pData(protocolData)) <- sns

 	protocolData(cnSet) <- protocolData

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

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

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

 	cnSet$SKW <- SKW

 	cnSet$SNR <- SNR

 	cnSet$gender <- genderff

@@ -192,7 +190,9 @@ snprmaAffy <- function(cnSet, filenames,

 	pkgname <- getCrlmmAnnotationName(annotation(cnSet))

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

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

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

+	sampleBatches <- splitIndicesByLength(seq_along(filenames), ocSamples()/getDoParWorkers())

+##	if(length(sampleBatches) < getDoParWorkers())

+##		sampleBatches <- splitIndicesByNode(seq_along(filenames))

 	ocLapply(sampleBatches,

 		 rsprocessCEL,

 		 filenames=filenames,

@@ -225,19 +225,25 @@ genotype <- function(filenames,

 		     recallRegMin=1000,

 		     gender=NULL,

 		     returnParams=TRUE,

-		     badSNP=0.7){

-	stopifnot(isPackageLoaded("ff"))

+		     badSNP=0.7,

+		     genome=c("hg19", "hg18")){

+	if(!isPackageLoaded("ff")) stop("ff package must be loaded")

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

+	if (any(duplicated(sns))) stop("sample identifiers must be unique")

+	genome <- match.arg(genome)

 	cnSet <- constructAffy(filenames=filenames,

+			       sns=sns,

 			       cdfName=cdfName,

-			       batch=batch, verbose=verbose)

+			       batch=batch, verbose=verbose,

+			       genome=genome)

 	if(!(is(A(cnSet), "ff") || is(A(cnSet), "ffdf"))){

 		stop("The ff package is required for this function.")

 	}

 	cnSet@mixtureParams <- snprmaAffy(cnSet, filenames=filenames,

-				    mixtureSampleSize=mixtureSampleSize,

-				    eps=eps,

-				    seed=seed,

-				    verbose=verbose)

+					  mixtureSampleSize=mixtureSampleSize,

+					  eps=eps,

+					  seed=seed,

+					  verbose=verbose)

 	ok <- cnrmaAffy(cnSet=cnSet, filenames=filenames,

 			cdfName=annotation(cnSet), seed=seed,

 			verbose=verbose)

@@ -697,6 +703,10 @@ fit.lm2 <- function(strata,

 	flags <- as.matrix(flags(object)[ii, batch.index])

 	fns <- featureNames(object)[ii]

 	fns.noflags <- fns[rowSums(flags, na.rm=T) == 0]

+	if(length(fns.noflags) == 0) {

+		warning("All features in one of the probeset batches were flagged. The data is processed in probeset batches to reduce RAM, but the error suggests that increasing the size of the batches may help.  For example, ocProbesets(4*ocProbesets()) would increase the current size of the probeset batches by 4.")

+		fns.noflags <- featureNames(object)[marker.index] ## allow processing to continue

+	}

 	snp.index <- sample(match(fns.noflags, featureNames(object)), 5000)

 	##

 	nuA.np <- as.matrix(nuA(object)[marker.index, batch.index])

@@ -1089,6 +1099,15 @@ fit.lm3 <- function(strata,

 	phiB(object)[marker.index, batch.index] <- phiB

 	phiPrimeA(object)[marker.index, batch.index] <- phiA2

 	phiPrimeB(object)[marker.index, batch.index] <- phiB2

+	##

+	if(enough.women){

+		medianA.AA(object)[marker.index, batch.index] <- do.call("cbind", lapply(medianA.Flist, function(x) x[, 1]))

+		medianA.AB(object)[marker.index, batch.index] <- do.call("cbind", lapply(medianA.Flist, function(x) x[, 2]))

+		medianA.BB(object)[marker.index, batch.index] <- do.call("cbind", lapply(medianA.Flist, function(x) x[, 3]))

+		medianB.AA(object)[marker.index, batch.index] <- do.call("cbind", lapply(medianB.Flist, function(x) x[, 1]))

+		medianB.AB(object)[marker.index, batch.index] <- do.call("cbind", lapply(medianB.Flist, function(x) x[, 2]))

+		medianB.BB(object)[marker.index, batch.index] <- do.call("cbind", lapply(medianB.Flist, function(x) x[, 3]))

+	}

 	if(is.lds) {close(object); return(TRUE)} else return(object)

 }

@@ -1145,6 +1164,10 @@ fit.lm4 <- function(strata,

 	fns <- featureNames(object)[ii]

 	flags <- as.matrix(flags(object)[ii, batch.index])

 	fns.noflags <- fns[rowSums(flags, na.rm=T) == 0]

+	if(length(fns.noflags) == 0){

+		warning("All features in one of the probeset batches were flagged. The data is processed in probeset batches to reduce RAM, but the error suggests that increasing the size of the batches may help.  For example, ocProbesets(4*ocProbesets()) would increase the current size of the probeset batches by 4.")

+		fns.noflags <- featureNames(object)[marker.index] ## allow processing to continue

+	}

 	snp.index <- sample(match(fns.noflags, featureNames(object)), 10000)

 	##

 	N.AA <- as.matrix(N.AA(object)[snp.index, batch.index, drop=FALSE])

@@ -1260,7 +1283,7 @@ cnrma2 <- function(A, filenames, row.names, verbose=TRUE, seed=1, cdfName, sns){

 	pkgname <- getCrlmmAnnotationName(cdfName)

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

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

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

+	sampleBatches <- splitIndicesByLength(seq_along(filenames), ocSamples()/getDoParWorkers())

 	if(verbose) message("Processing nonpolymorphic probes for ", length(filenames), " files.")

 	## updates A

 	ocLapply(sampleBatches,

@@ -1540,10 +1563,10 @@ shrinkSummary <- function(object,

 		marker.index <- whichMarkers(type[[1]], is.snp, is.autosome, is.annotated, is.X)

 	}

 	if(is.lds){

-		##index.list <- splitIndicesByLength(marker.index, ocProbesets())

-		if(parStatus()){

-			index.list <- splitIndicesByNode(marker.index)

-		} else index.list <- splitIndicesByLength(marker.index, ocProbesets())

+		index.list <- splitIndicesByLength(marker.index, ocProbesets())

+##		if(parStatus()){

+##			index.list <- splitIndicesByNode(marker.index)

+##		} else index.list <- splitIndicesByLength(marker.index, ocProbesets())

 		ocLapply(seq(along=index.list),

 			 shrinkGenotypeSummaries,

 			 index.list=index.list,

@@ -1594,10 +1617,10 @@ genotypeSummary <- function(object,

 	myf <- summaryFxn(type[[1]])

 	FUN <- get(myf)

 	if(is.lds){

-		##index.list <- splitIndicesByLength(marker.index, ocProbesets())

-		if(parStatus()){

-			index.list <- splitIndicesByNode(marker.index)

-		} else index.list <- splitIndicesByLength(marker.index, ocProbesets())

+		index.list <- splitIndicesByLength(marker.index, ocProbesets())

+##		if(parStatus()){

+##			index.list <- splitIndicesByNode(marker.index)

+##		} else index.list <- splitIndicesByLength(marker.index, ocProbesets())

 		ocLapply(seq(along=index.list),

 			 FUN,

 			 index.list=index.list,

@@ -2116,9 +2139,10 @@ estimateCnParameters <- function(object,

 	myfun <- lmFxn(type[[1]])

 	FUN <- get(myfun)

 	if(is.lds){

-		if(parStatus()){

-			index.list <- splitIndicesByNode(marker.index)

-		} else index.list <- splitIndicesByLength(marker.index, ocProbesets())

+		index.list <- splitIndicesByLength(marker.index, ocProbesets())

+##		if(parStatus()){

+##			index.list <- splitIndicesByNode(marker.index)

+##		} else index.list <- splitIndicesByLength(marker.index, ocProbesets())

 		ocLapply(seq(along=index.list),

 			 FUN,

 			 index.list=index.list,

@@ -2784,89 +2808,73 @@ ABpanel <- function(x, y, predictRegion,

 	}

 }

-calculateRTheta <- function(cnSet, genotype=c("AA", "AB", "BB"), batch.name){

-	genotype <- match.arg(genotype)

-	j <- match(batch.name, batchNames(cnSet))

-	##centers <- medians(cnSet, i=snp.index, j)

-	centers <- medians(cnSet, i=seq_len(nrow(cnSet)), j)

-	theta <- matrix(NA, nrow(cnSet), 2)

-	colnames(theta) <- c("theta", "R")

-	x <- centers[, "A", genotype, j]

-	y <- centers[, "B", genotype, j]

-	## small imputed values -- should fix imputeCenter

-	x[x < 64] <- 64

-	y[y < 64] <- 64

-	theta[, "theta"] <- atan2(y, x)*2/pi

-##	if(any(is.na(theta[, "theta"]))){

-##		##stop("NA's in thetas.  Can not compute theta / R values for batches with fewer than 10 samples")

-##	}

-	## so that 'R' is available for NP probes

-	y[is.na(y)] <- 0

-	theta[, "R"] <- x+y

-	theta[!isSnp(cnSet), 1] <- 1

-	if(any(is.na(theta[, "theta"]))){

-		stop("NA's in thetas.  Can not compute theta / R values for batches with fewer than 10 samples")

-	}

-	return(theta)

-}

-

-

-

-

-calculateBAFandLRR <- function(cnSet) {

-  ## Calculates B allele frequency and log2 R ratio for all snps for a given cnSet

-  ## Returns a 3-D array of size Num Of Snps x Num of Samples x 2

-  ## where result[,,1] is B allele frequency and result[,,2] is log2 R ratio

-

-  snp.index <- which(isSnp(cnSet))

-  b <- (B(cnSet))[snp.index, ]

-  a <- (A(cnSet))[snp.index, ]

-  centers<-medians(cnSet, i=snp.index, j=1:length(unique(batch(cnSet))))

-  bafAndLrr <- array(NA_integer_, dim=c(length(snp.index), length(sampleNames(cnSet)), 2), dimnames=list(featureNames(cnSet)[snp.index], sampleNames(cnSet),c("baf", "lrr")))

-

-  for (batch in unique(batch(cnSet))) {

-

-    ## get batch- and snp-specific centers

-    center.aa.x <- centers[, 1, 1, batch]

-    center.aa.y <- centers[, 2, 1, batch]

-    center.ab.x <- centers[, 1, 2, batch]

-    center.ab.y <- centers[, 2, 2, batch]

-    center.bb.x <- centers[, 1, 3, batch]

-    center.bb.y <- centers[, 2, 3, batch]

-    theta.aa <- atan2(center.aa.y, center.aa.x) * 2 / pi

-    r.aa <- center.aa.x + center.aa.y

-    theta.ab <- atan2(center.ab.y, center.ab.x) * 2 / pi

-    r.ab <- center.ab.x + center.ab.y

-    theta.bb <- atan2(center.bb.y, center.bb.x) * 2 / pi

-    r.bb <- center.bb.x + center.bb.y

-

-    index.sample <- which(batch(cnSet) %in% batch)

-    for (i in index.sample) {

-      theta <- atan2(b[, i], a[, i]) * 2 / pi

-      r <- a[, i] + b[, i]

-

-      baf <- vector(mode="numeric", length=length(theta))

-      baf[theta < theta.aa] <- 0

-      baf[theta >= theta.bb] <- 1

-      idx.ab <- which(theta >= theta.aa & theta < theta.ab)

-      baf[idx.ab] <- .5 * (theta[idx.ab] - theta.aa[idx.ab]) / (theta.ab[idx.ab] - theta.aa[idx.ab])

-      idx.bb <- which(theta >= theta.ab & theta < theta.bb)

-      baf[idx.bb] <- .5 * (theta[idx.bb] - theta.ab[idx.bb]) / (theta.bb[idx.bb] - theta.ab[idx.bb]) + .5

-      bafAndLrr[, i ,1] <- baf

-

-      r.expected <- vector(mode="numeric", length=length(r))

-      r.expected[theta < theta.aa] <- r.aa[theta < theta.aa]

-      r.expected[theta >= theta.bb] <- r.bb[theta >= theta.bb]

-

-      idx.ab <- which(theta >= theta.aa & theta < theta.ab)

-      r.expected[idx.ab] <- (theta[idx.ab] - theta.aa[idx.ab]) * (r.ab[idx.ab]-r.aa[idx.ab]) /

-        (theta.ab[idx.ab] - theta.aa[idx.ab]) + r.aa[idx.ab]

-      idx.bb <- which(theta >= theta.ab & theta < theta.bb)

-      r.expected[idx.bb] <- (theta[idx.bb] - theta.ab[idx.bb]) * (r.bb[idx.bb]-r.ab[idx.bb])/

-        (theta.bb[idx.bb] - theta.ab[idx.bb]) + r.ab[idx.bb]

-      bafAndLrr[, i, 2] <- log2(r / r.expected)

-

-    }

-  }

-  return(bafAndLrr)

+calculateRTheta <- function(object, ##genotype=c("AA", "AB", "BB"),

+			    batch.name,

+			    feature.index){

+	index.np <- which(!(isSnp(object)[feature.index]))

+	j <- match(batch.name, batchNames(object))

+	if(missing(j)) stop("batch.name did not match in batchNames(object)")

+	CHR <- unique(chromosome(object)[feature.index])

+	isX <- CHR == "X"

+	centers <- getMedians(batchStatistics(object))

+	lapply(centers, open)

+	centersMatrix <- lapply(centers, function(x, i, j) x[i, j], j=j, i=feature.index)

+	lapply(centers, close)

+	centersMatrix <- do.call(cbind, centersMatrix)

+	centersA <- centersMatrix[, 1:3]

+	centersB <- centersMatrix[, 4:6]

+	rm(centers, centersMatrix); gc()

+	centersA[centersA < 64] <- 64

+	centersB[centersB < 64] <- 64

+	theta <- atan2(centersB, centersA) * 2/pi

+	centersA[is.na(centersA)] <- 0

+	centersB[is.na(centersB)] <- 0

+	if(length(index.np) > 0) theta[index.np, ] <- 1

+	##

+	r <- centersA+centersB

+	J <- which(batch(object)==batch.name)

+	open(A(object))

+	open(B(object))

+	a <- as.matrix(A(object)[feature.index, J, drop=FALSE])

+	b <- as.matrix(B(object)[feature.index, J, drop=FALSE])

+	close(A(object))

+	close(B(object))

+	if(length(index.np) > 0) b[index.np, ] <- 0L

+	obs.theta <- atan2(b, a)*2/pi

+	calculateBandR <- function(o, r, theta, not.na, M){

+		lessAA <- o < theta[, 1]

+		lessAB <- o < theta[, 2]

+		lessBB <- o < theta[, 3]

+		ii <- which(lessAA & not.na)

+		jj <- which(!lessBB & not.na)

+		i <- which(!lessAA & lessAB & not.na)

+		j <- which(!lessAB & lessBB & not.na)

+		M[i, 1] <- 0.5*(o[i]-theta[i,1])/(theta[i,2]-theta[i,1])

+		M[j, 1] <- 0.5*(o[j]-theta[j,2])/(theta[j,3]-theta[j,2]) + 0.5

+		M[ii, 1] <- 0

+		M[jj, 1] <- 1

+		M[i, 2] <- (o[i]-theta[i,1])*(r[i,2]-r[i,1])/(theta[i,2]-theta[i,1]) + r[i, 1]

+		M[j, 2] <- (o[j]-theta[j,2])*(r[j,3]-r[j,2])/(theta[j,3]-theta[j,2]) + r[j, 2]

+		M[ii, 2] <- r[ii, 1]

+		M[jj, 2] <- r[jj, 3]

+		return(M)

+	}

+	not.na <- !is.na(theta[,1])

+	r.expected <- bf <- matrix(NA, length(feature.index), ncol(a))

+	M <- matrix(NA, length(feature.index), 2)

+	for(j in seq_len(ncol(a))){

+		tmp <- calculateBandR(obs.theta[, j], r=r, theta=theta, not.na=not.na, M=M)

+		bf[, j] <- tmp[,1]

+		r.expected[,j] <- tmp[,2]

+	}

+	bf[bf < 0] <- 0

+	bf[bf > 1] <- 1

+	if(length(index.np) > 0) r.expected[index.np, ] <- centersA[index.np, 1]

+	obs.r <- a+b

+	lrr <- log2(obs.r/r.expected)

+	dimnames(bf) <- dimnames(lrr) <- list(featureNames(object)[feature.index],

+					      sampleNames(object)[J])

+	bf <- integerMatrix(bf, 1000)

+	lrr <- integerMatrix(lrr, 100)

+	return(list(baf=bf, lrr=lrr))

 }

@@ -7,6 +7,7 @@ crlmm <- function(filenames, row.names=TRUE, col.names=TRUE,

   if ((load.it | save.it) & missing(intensityFile))

 	  stop("'intensityFile' is missing, and you chose either load.it or save.it")

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

+  if (any(duplicated(sns))) stop("sample identifiers are not unique")

   if (!missing(intensityFile))

     if (load.it & !file.exists(intensityFile)){

       load.it <- FALSE

@@ -312,24 +313,26 @@ crlmm2 <- function(filenames, row.names=TRUE, col.names=TRUE,

 }

 imputeGender <- function(A, B, XIndex, YIndex, SNR, SNRMin){

-	##XMedian <- apply(log2(A[XIndex,,

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

-	a <- log2(A[XIndex,,drop=FALSE])

-	b <- log2(B[XIndex,,drop=FALSE])

-	meds.X <- (apply(a+b, 2, median))/2

-	##YIndex <- which(isSnp(gtSet) & chromosome(gtSet)==24)

-	a <- log2(A[YIndex,,drop=FALSE])

-	b <- log2(B[YIndex,,drop=FALSE])

-	meds.Y <- (apply(a+b, 2, median))/2

-	R <- meds.X - meds.Y

-	##SNR <- gtSet$SNR[]

-	##SNRMin=5

-	if(sum(SNR > SNRMin) == 1){

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

-		gender <- ifelse(R[SNR[] > SNRMin] > 0.5, 2L, 1L)

-	} else{

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

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

+	if(length(YIndex) > 0){

+		a <- log2(A[XIndex,,drop=FALSE])

+		b <- log2(B[XIndex,,drop=FALSE])

+		meds.X <- (apply(a+b, 2, median))/2

+		a <- log2(A[YIndex,,drop=FALSE])

+		b <- log2(B[YIndex,,drop=FALSE])

+		meds.Y <- (apply(a+b, 2, median))/2

+		R <- meds.X - meds.Y

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

+			gender <- ifelse(R[SNR[] > SNRMin] > 0.5, 2L, 1L)

+		} else{

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

+		}

+	} else {

+		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"]]

+		}

 	}

 	return(gender)

 }

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

-# function below works OK provided all .idat files are in the current working directory

+					# function below works OK provided all .idat files are in the current working directory

 # - could add an option to allow files in Illumina directory structure to be handled

 # or to use the optional 'Path' column in sampleSheet

 # - there is a lot of header information that is currently discarded - could try and store this somewhere in the resulting NChannelSet

@@ -36,7 +36,7 @@ readIdatFiles = function(sampleSheet=NULL,

 		       }

 	       }

 	       pd = new("AnnotatedDataFrame", data = sampleSheet)

-	       sampleNames(pd) = basename(arrayNames)

+	       sampleNames(pd) = make.unique(basename(arrayNames))

        }

        if(is.null(arrayNames)) {

                arrayNames = gsub(paste(sep, fileExt$green, sep=""), "", dir(pattern=fileExt$green, path=path))

@@ -95,7 +95,7 @@ readIdatFiles = function(sampleSheet=NULL,

 	       if(i==1) {

 		       if(is.null(ids) && !is.null(G)){

 			       ids = idsG

-		       } else stop("Could not find probe IDs")

+		       } # else stop("Could not find probe IDs")

 		       nprobes = length(ids)

 		       narrays = length(arrayNames)

 		       RG = new("NChannelSet",

@@ -106,8 +106,8 @@ readIdatFiles = function(sampleSheet=NULL,

 				 phenoData=pd, storage.mode="environment")

 		       featureNames(RG) = ids

 		       if(!is.null(sampleSheet) && !is.null(sampleSheet$Sample_ID)){

-		            sampleNames(RG) = sampleSheet$Sample_ID

-		       } else  sampleNames(RG) = arrayNames

+		            sampleNames(RG) = make.unique(sampleSheet$Sample_ID)

+		       } else  sampleNames(RG) = make.unique(basename(arrayNames))

 		       gc(verbose=FALSE)

 	       }

 	       if(length(ids)==length(idsG)) {

@@ -524,7 +524,6 @@ readGenCallOutput = function(file, path=".", cdfName,

 RGtoXY = function(RG, chipType, verbose=TRUE) {

-

   needToLoad <- !all(sapply(c('addressA', 'addressB', 'base'), isLoaded))

   if(needToLoad){

 	  chipList = c("human1mv1c",# 1M

@@ -537,12 +536,14 @@ RGtoXY = function(RG, chipType, verbose=TRUE) {

 	  "human1mduov3b",          # 1M Duo

 	  "humanomni1quadv1b",      # Omni1 quad

 	  "humanomni25quadv1b",     # Omni2.5 quad

+	  "humanomni258v1a",        # Omni2.5 8 sample

 	  "humanomniexpress12v1b",  # Omni express 12

 	  "humanimmuno12v1b",       # Immuno chip 12

           "humancytosnp12v2p1h")    # CytoSNP 12

+	  ## RS: added cleancdfname()

 	  if(missing(chipType)){

-		  chipType = match.arg(annotation(RG), chipList)

-	  } else chipType = match.arg(chipType, chipList)

+		  chipType = match.arg(cleancdfname(annotation(RG)), chipList)

+	  } else chipType = match.arg(cleancdfname(chipType), chipList)

 	  pkgname = getCrlmmAnnotationName(chipType)

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

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

@@ -1106,7 +1107,7 @@ getProtocolData.Illumina = function(filenames, sep="_", fileExt="Grn.idat", verb

                          Position = rep(NA, narrays))

        scanDates = data.frame(ScanDate=rep(NA, narrays), DecodeDate=rep(NA, narrays))

-       rownames(scanDates) = gsub(paste(sep, fileExt, sep=""), "", filenames)

+       rownames(scanDates) <- make.unique(gsub(paste(sep, fileExt, sep=""), "", filenames))

        ## read in the data

        for(i in seq_along(filenames)) {

                if(verbose)

@@ -1132,7 +1133,7 @@ getProtocolData.Illumina = function(filenames, sep="_", fileExt="Grn.idat", verb

 			    data=scanDates,

 			    varMetadata=data.frame(labelDescription=colnames(scanDates),

 			                           row.names=colnames(scanDates)))

-	return(protocoldata)

+       return(protocoldata)

 }

@@ -1171,7 +1172,7 @@ constructInf <- function(sampleSheet=NULL,

 			}

 		}

 		pd = new("AnnotatedDataFrame", data = sampleSheet)

-		sampleNames(pd) = basename(arrayNames)

+		sampleNames(pd) = make.unique(basename(arrayNames))

 	}

 	if(is.null(arrayNames)) {

 		arrayNames = gsub(paste(sep, fileExt$green, sep=""), "", dir(pattern=fileExt$green, path=path))

@@ -1210,17 +1211,29 @@ constructInf <- function(sampleSheet=NULL,

 	if(verbose) message("Initializing container for genotyping and copy number estimation")

 	featureData = getFeatureData(cdfName, copynumber=TRUE)

 	nr = nrow(featureData); nc = narrays

+	sns <-  if (!is.null(sampleSheet) && !is.null(sampleSheet$Sample_ID)) {

+		make.unique(sampleSheet$Sample_ID)

+        } else{

+		make.unique(basename(arrayNames))

+	}

+	biga <- initializeBigMatrix(name="A", nr, nc)

+	bigb <- initializeBigMatrix(name="B", nr, nc)

+	bigc <- initializeBigMatrix(name="call", nr, nc)

+	bigd <- initializeBigMatrix(name="callPr", nr,nc)

+	colnames(biga) <- colnames(bigb) <- colnames(bigc) <- colnames(bigd) <- sns

 	cnSet <- new("CNSet",

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

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

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

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

+		     alleleA=biga,

+		     alleleB=bigb,

+		     call=bigc,

+		     callProbability=bigd,

 		     annotation=cdfName,

 		     featureData=featureData,

 		     batch=batch)

-        if (!is.null(sampleSheet) && !is.null(sampleSheet$Sample_ID)) {

-	   	sampleNames(cnSet) = sampleSheet$Sample_ID

-        } else  sampleNames(cnSet) = arrayNames

+##        if (!is.null(sampleSheet) && !is.null(sampleSheet$Sample_ID)) {

+##		sampleNames(cnSet) = make.unique(sampleSheet$Sample_ID)

+##        } else{

+##		sampleNames(cnSet) <- make.unique(basename(arrayNames))

+##	}

 	if(saveDate){

 		protocolData = getProtocolData.Illumina(grnidats, sep=sep, fileExt=fileExt$green, verbose=verbose)

 	} else{

@@ -1233,7 +1246,7 @@ constructInf <- function(sampleSheet=NULL,

 	cnSet$gender <- initializeBigVector("gender", ncol(cnSet), vmode="integer")

 	cnSet$SNR = initializeBigVector("crlmmSNR-", ncol(cnSet), "double")

 	cnSet$SKW = initializeBigVector("crlmmSKW-", ncol(cnSet), "double")

-	sampleNames(cnSet) <- basename(sampleNames(cnSet))

+	##sampleNames(cnSet) <- basename(sampleNames(cnSet))

 	return(cnSet)

 }

 construct.Illumina <- constructInf

@@ -1263,8 +1276,7 @@ preprocessInf <- function(cnSet,

 	is.snp = isSnp(cnSet)

 	snp.index = which(is.snp)

 	narrays = ncol(cnSet)

-#	sampleBatches <- splitIndicesByLength(seq(length=ncol(cnSet)), ocSamples())

-	sampleBatches <- splitIndicesByNode(seq(length=ncol(cnSet)))

+	sampleBatches <- splitIndicesByLength(seq_len(ncol(cnSet)), ocSamples()/getDoParWorkers())

 	mixtureParams = initializeBigMatrix("crlmmMixt-", 4, narrays, "double")

 	ocLapply(seq_along(sampleBatches),

 		 crlmm:::processIDAT,

@@ -50,6 +50,21 @@ setMethod("medians", signature(object="AssayData"),

 		  return(res)

 })

+getMedians <- function(object){

+	medianA.AA <- assayDataElement(object, "medianA.AA")

+	medianA.AB <- assayDataElement(object, "medianA.AB")

+	medianA.BB <- assayDataElement(object, "medianA.BB")

+	medianB.AA <- assayDataElement(object, "medianB.AA")

+	medianB.AB <- assayDataElement(object, "medianB.AB")

+	medianB.BB <- assayDataElement(object, "medianB.BB")

+	list(A.AA=medianA.AA,

+	     A.AB=medianA.AB,

+	     A.BB=medianA.BB,

+	     B.AA=medianB.AA,

+	     B.AB=medianB.AB,

+	     B.BB=medianB.BB)

+}

+

 setMethod("mads", signature(object="AssayData"),

 	  function(object, ...){

 		  batchnames <- sampleNames(object)

@@ -1,63 +1,6 @@

 setMethod("posteriorMean", signature(object="CNSet"), function(object) assayDataElement(object, "posteriorMean"))

 setReplaceMethod("posteriorMean", signature(object="CNSet", value="matrix"), function(object, value) assayDataElementReplace(object, "posteriorMean", value))

-setAs("CNSet", "oligoSnpSet", function(from, to){

-	cnSet2oligoSnpSet(from)

-})

-

-cnSet2oligoSnpSet <- function(object){

-	row.index <- seq_len(nrow(object))

-	col.index <- seq_len(ncol(object))

-	is.lds <- ifelse(is(calls(object), "ff_matrix") | is(calls(object), "ffdf"), TRUE, FALSE)

-	if(is.lds) stopifnot(isPackageLoaded("ff"))

-	b.r <- calculateRBaf(object)

-	b <- integerMatrix(b.r[[1]], 1000)

-	r <- integerMatrix(b.r[[2]], 100)

-##	if(is.lds){

-##		## initialize a big matrix for raw copy number

-##		message("creating an ff object for storing total copy number")

-##		tcn <- initializeBigMatrix(name="total_cn", nrow(object), ncol(object), vmode="double")

-##		for(j in 1:ncol(object)){

-##			tcn[, j] <- totalCopynumber(object, i=row.index, j=j)

-##		}

-##	} else {

-##		if(ncol(object) > 5){

-##			##this can be memory intensive, so we try to be careful

-##			col.index <- splitIndicesByLength(seq(length=ncol(object)), 5)

-##			tcn <- matrix(NA, nrow(object), ncol(object))

-##			dimnames(tcn) <- list(featureNames(object), sampleNames(object))

-##			rows <- 1:nrow(object)

-##			for(i in seq_along(col.index)){

-##				cat(".")

-##				j <- col.index[[i]]

-##				cnSet <- object[, j]

-##				tcn[, j] <- totalCopynumber(cnSet, i=row.index, j=1:ncol(cnSet))

-##				rm(cnSet); gc()

-##			}

-##			cat("\n")

-##		} else {

-##			tcn <- totalCopynumber(object, i=row.index, j=col.index)

-##		}

-##	}

-##	message("Transforming copy number to log2 scale")

-##	tcn[tcn < 0.1] <- 0.1

-##	tcn[tcn > 8] <- 8

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

-	tmp <- new("oligoSnpSet",

-		   copyNumber=integerMatrix(b.r[[2]], scale=100),

-		   call=calls(object),

-		   callProbability=snpCallProbability(object),

-		   annotation=annotation(object),

-		   featureData=featureData(object),

-		   phenoData=phenoData(object),

-		   experimentData=experimentData(object),

-		   protocolData=protocolData(object),

-		   is.integer=FALSE)

-	tmp <- assayDataElementReplace(tmp, "baf", integerMatrix(b.r[[1]], scale=1000))

-	return(tmp)

-}

-

-

 linearParamElementReplace <- function(obj, elt, value) {

     storage.mode <- storageMode(batchStatistics(obj))

     switch(storage.mode,

@@ -185,7 +128,7 @@ ACN <- function(object, allele, i , j){

 	}

 	## Define batch.index and sample.index

 	if(!missing.j) {

-		batches <- unique(as.character(batch(object))[j])

+		batches <- unique(as.character(batch(object)[j]))

 		##batches <- as.character(batch(object)[j])

 		batch.index <- match(batches, batchNames(object))

 	} else {

@@ -587,78 +530,129 @@ setMethod("xyplot", signature(x="formula", data="CNSet"),

 })

 setMethod("calculateRBaf", signature(object="CNSet"),

-	  function(object, batch.name){

-		  all.autosomes <- all(chromosome(object) < 23)

-		  if(!all.autosomes){

-			  stop("method currently only defined for chromosomes 1-22")

-		  }

-		  if(missing(batch.name)) batch.name <- batchNames(object)[1]

-		  stopifnot(batch.name %in% batchNames(object))

-		  if(length(batch.name) > 1){

-			  warning("only the first batch in batch.name processed")

-			  batch.name <- batch.name[1]

-		  }

-		  RTheta.aa <- calculateRTheta(object, "AA", batch.name)

-		  RTheta.ab <- calculateRTheta(object, "AB", batch.name)

-		  RTheta.bb <- calculateRTheta(object, "BB", batch.name)

-

-		  J <- which(batch(object) == batch.name)

-

-		  theta.aa <- matrix(RTheta.aa[, "theta"], nrow(object), length(J), byrow=FALSE)

-		  theta.ab <- matrix(RTheta.ab[, "theta"], nrow(object), length(J), byrow=FALSE)

-		  theta.bb <- matrix(RTheta.bb[, "theta"], nrow(object), length(J), byrow=FALSE)

-

-		  a <- A(object)[, J, drop=FALSE]

-		  b <- B(object)[, J, drop=FALSE] ## NA's for b where nonpolymorphic

-		  is.np <- !isSnp(object)

-		  ##b[is.np, ] <- a[is.np, ]

-		  b[is.np, ] <- 0L

-		  dns <- dimnames(a)

-		  dimnames(a) <- dimnames(b) <- NULL

-		  obs.theta <- atan2(b, a)*2/pi

-

-		  lessAA <- obs.theta < theta.aa

-		  lessAB <- obs.theta < theta.ab

-		  lessBB <- obs.theta < theta.bb

-		  grAA <- !lessAA

-		  grAB <- !lessAB

-		  grBB <- !lessBB

-		  not.na <- !is.na(theta.aa)

-		  I1 <- grAA & lessAB & not.na

-		  I2 <- grAB & lessBB & not.na

-

-		  bf <- matrix(NA, nrow(object), ncol(a))

-		  bf[I1] <- 0.5 * ((obs.theta-theta.aa)/(theta.ab-theta.aa))[I1]

-		  bf[I2] <- (.5 * (obs.theta - theta.ab) / (theta.bb - theta.ab))[I2] + 0.5

-		  bf[lessAA] <- 0

-		  bf[grBB] <- 1

-

-		  r.expected <- matrix(NA, nrow(object), ncol(a))

-		  r.aa <- matrix(RTheta.aa[, "R"], nrow(object), length(J), byrow=FALSE)

-		  r.ab <- matrix(RTheta.ab[, "R"], nrow(object), length(J), byrow=FALSE)

-		  r.bb <- matrix(RTheta.bb[, "R"], nrow(object), length(J), byrow=FALSE)

-		  rm(RTheta.aa, RTheta.ab, RTheta.bb); gc()

-		  obs.r <- a+b

-

-		  lessAA <- lessAA & not.na

-		  grBB <- grBB & not.na

-		  tmp <- ((obs.theta - theta.aa) * (r.ab-r.aa)/(theta.ab-theta.aa))[I1] + r.aa[I1]

-		  r.expected[I1] <- tmp

-		  tmp <- ((obs.theta - theta.ab) * (r.bb - r.ab)/(theta.bb-theta.ab))[I2] + r.ab[I2]

-		  r.expected[I2] <- tmp

-		  r.expected[lessAA] <- r.aa[lessAA]

-		  r.expected[grBB] <- r.bb[grBB]

+	  function(object, batch.name, chrom){

+		  calculateRBafCNSet(object, batch.name, chrom)

+	  })

-		  index.np <- which(is.np)

-		  if(length(index.np) > 0){

-			  a.np <- A(object)[index.np, J]

-			  meds <- rowMedians(a.np, na.rm=TRUE)

-			  r.expected[index.np, ] <- matrix(meds, length(index.np), ncol(a.np))

-		  }

-		  lrr <- log2(obs.r/r.expected)

+calculateRBafCNSet <- function(object, batch.name, chrom){

+	if(missing(batch.name)) batch.name <- batchNames(object)

+	if(missing(chrom)) chrom <- unique(chromosome(object))

+	if(!(all(batch.name %in% batchNames(object)))) stop("batch.name must be belong to batchNames(object)")

+	chr <- chromosome(object)

+	valid.chrs <- chr <= 23 & chr %in% chrom

+	index <- which(valid.chrs)

+	indexlist <- split(index, chr[index])

+	J <- which(batch(object) %in% batch.name)

+	sns <- sampleNames(object)[J]

+	sampleindex <- split(J, batch(object)[J])

+	if(!all(valid.chrs)) warning("Only computing log R ratios and BAFs for autosomes and chr X")

+	## if ff package is loaded, these will be ff objects

+	chr <- names(indexlist)

+	rlist <- blist <- vector("list", length(indexlist))

+	for(i in seq_along(indexlist)){

+		I <- indexlist[[i]]

+		nr <- length(I)

+		CHR <- names(indexlist)[i]

+		bafname <- paste("baf_chr", CHR, sep="")

+		rname <- paste("lrr_chr", CHR, sep="")

+		bmatrix <- initializeBigMatrix(bafname, nr=nr, nc=length(sns), vmode="integer")

+		rmatrix <- initializeBigMatrix(rname, nr=nr, nc=length(sns), vmode="integer")

+		colnames(rmatrix) <- colnames(bmatrix) <- sns

+		## put rownames in order of physical position

+		ix <- order(position(object)[I])

+		I <- I[ix]

+		rownames(rmatrix) <- rownames(bmatrix) <- featureNames(object)[I]

+		for(j in seq_along(sampleindex)){

+			bname <- batch.name[j]

+			J <- sampleindex[[j]]

+			res <- calculateRTheta(object=object,

+					       batch.name=bname,

+					       feature.index=I)

+			k <- match(sampleNames(object)[J], sns)

+			bmatrix[, k] <- res[["baf"]]

+			rmatrix[, k] <- res[["lrr"]]

+			##colnames(bmatrix)[J] <- colnames(rmatrix)[J] <- sampleNames(object)[J]

+		}

+		blist[[i]] <- bmatrix

+		rlist[[i]] <- rmatrix

+	}

+	res <- list(baf=blist,

+		    lrr=rlist)

+	return(res)

+}

-		  dimnames(bf) <- dimnames(lrr) <- dns

-		  res <- list(baf=bf,

-			      lrr=lrr)

-		  return(res)

-	  })

+##setMethod("calculateRBaf", signature(object="CNSet"),

+##	  function(object, batch.name){

+##		  all.autosomes <- all(chromosome(object) < 23)

+##		  if(!all.autosomes){

+##			  stop("method currently only defined for chromosomes 1-22")

+##		  }

+##		  if(missing(batch.name)) batch.name <- batchNames(object)[1]

+##		  stopifnot(batch.name %in% batchNames(object))

+##		  if(length(batch.name) > 1){

+##			  warning("only the first batch in batch.name processed")

+##			  batch.name <- batch.name[1]

+##		  }

+##		  RTheta.aa <- calculateRTheta(object, "AA", batch.name)

+##		  RTheta.ab <- calculateRTheta(object, "AB", batch.name)

+##		  RTheta.bb <- calculateRTheta(object, "BB", batch.name)

+##

+##		  J <- which(batch(object) == batch.name)

+##

+##		  theta.aa <- matrix(RTheta.aa[, "theta"], nrow(object), length(J), byrow=FALSE)

+##		  theta.ab <- matrix(RTheta.ab[, "theta"], nrow(object), length(J), byrow=FALSE)

+##		  theta.bb <- matrix(RTheta.bb[, "theta"], nrow(object), length(J), byrow=FALSE)

+##

+##		  a <- A(object)[, J, drop=FALSE]

+##		  b <- B(object)[, J, drop=FALSE] ## NA's for b where nonpolymorphic