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

 2010-04-24 B Carvalho committed version 1.7.1

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

+2010-05-25 R. Scharpf committed version 1.7.2

+** added example dataset

+** fixed vignette: inst/scripts/copynumber.Rnw

@@ -1,8 +1,8 @@

 Package: crlmm

 Type: Package

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

-Version: 1.7.1

-Date: 2010-04-16

+Version: 1.7.2

+Date: 2010-05-25

 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>

 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

@@ -23,7 +23,8 @@ Imports: affyio (>= 1.15.2),

          utils

 Suggests: hapmapsnp6,

           genomewidesnp6Crlmm (>= 1.0.2),

-          snpMatrix

+          snpMatrix,

+	  ellipse

 Collate: AllGenerics.R

 	 AllClasses.R

 	 methods-CNSet.R

@@ -33,6 +34,7 @@ Collate: AllGenerics.R

          crlmm-functions.R

          crlmm-illumina.R

 	 snprma-functions.R

+	 plot-methods.R

          utils.R

          zzz.R

 LazyLoad: yes

@@ -70,6 +70,7 @@ export(crlmm,

        batch,

        crlmmCopynumber2)

 ##export(initializeParamObject, biasAdjNP2)

+##export(construct)

@@ -8,3 +8,16 @@ setMethod("initialize", "CNSetLM", function(.Object, lM=new("list"), ...){

 	.Object@lM <- lM

 	.Object <- callNextMethod(.Object, ...)

 })

+

+setValidity("CNSetLM",

+	    function(object){

+		    if(!"batch" %in% varLabels(protocolData(object)))

+			    return("'batch' not defined in protocolData")

+		    if(!"chromosome" %in% fvarLabels(object))

+			    return("'chromosome' not defined in featureData")

+		    if(!"position" %in% fvarLabels(object))

+			    return("'position' not defined in featureData")

+		    if(!"isSnp" %in% fvarLabels(object))

+			    return("'isSnp' not defined in featureData")

+		    return(TRUE)

+	    })

@@ -24,7 +24,6 @@ getFeatureData.Affy <- function(cdfName, copynumber=FALSE){

 	gns <- getVarInEnv("gns")

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

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

-

 	if(copynumber){

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

 		cnProbes <- get("cnProbes")

@@ -40,7 +39,6 @@ getFeatureData.Affy <- function(cdfName, copynumber=FALSE){

 	M[index, "isSnp"] <- 1L

 	index <- which(is.na(M[, "isSnp"]))

 	M[index, "isSnp"] <- 1L

-

 	if(copynumber){

 		index <- match(rownames(cnProbes), rownames(M)) #only snp probes in M get assigned position

 		M[index, "position"] <- cnProbes[, grep("pos", colnames(cnProbes))]

@@ -54,29 +52,42 @@ getFeatureData.Affy <- function(cdfName, copynumber=FALSE){

 	##crlmmOpts$snpRange <- range(snpIndex)

 	##crlmmOpts$npRange <- range(npIndex)

 }

+

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

-		      sns, verbose=TRUE){

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

+		      sns, verbose=TRUE, batch){

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

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

 	protocolData <- getProtocolData.Affy(filenames)

+	if(missing(batch)){

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

+	} else 	{

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

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

+		protocolData$batch <- batch

+	}

+	rownames(pData(protocolData)) <- NULL

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

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

-	callSet <- new("SnpSuperSet", 

+	callSet <- new("CNSetLM", 

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

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

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

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

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

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

+		       protocolData=protocolData,

 		       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)

-	protocolData(callSet) <- protocolData

-	sampleNames(callSet) <- sns

+	rownames(pData(callSet)) <- NULL

+	phenoData(callSet)$sampleNames <- sns

 	return(callSet)

 }

+

 genotype <- function(filenames,

 		     cdfName,

 		     batch,

@@ -97,27 +108,20 @@ genotype <- function(filenames,

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

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

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

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

-	}	

 	## callSet contains potentially very big matrices

 	## More big matrices are created within snprma, that will then be removed.

 	callSet <- construct(filenames=filenames,

 			     cdfName=cdfName,

 			     copynumber=copynumber,

 			     sns=sns,

-			     verbose=verbose)

-	if(missing(batch)){

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

-	}

+			     verbose=verbose,

+			     batch=batch)

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

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

 	batches <- splitIndicesByLength(1:ncol(callSet), ocSamples())

 	iter <- 1

-	for(j in batches){

+##	for(j in batches){

+	j <- unlist(batches)

 		if(verbose) message("Batch ", iter, " of ", length(batches))

 		snprmaRes <- snprma(filenames=filenames[j],

 				    mixtureSampleSize=mixtureSampleSize,

@@ -167,7 +171,7 @@ genotype <- function(filenames,

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

 		callSet$gender[j] <- tmp$gender

 		iter <- iter+1

-	}

+##	}

 	return(callSet)

 }

@@ -208,23 +212,15 @@ genotype2 <- function(filenames,

 	}

 	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

+			     verbose=verbose,

+			     batch=batch)

+	##lM(callSet) <- initializeParamObject(list(featureNames(callSet), unique(protocolData(callSet)$batch)))

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

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

 	snprmaRes <- snprma2(filenames=filenames,

@@ -242,16 +238,23 @@ genotype2 <- function(filenames,

 	open(snprmaRes[["SKW"]])

 	open(snprmaRes[["mixtureParams"]])

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

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

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

+	##batches <- splitIndicesByLength(1:nrow(snprmaRes[["A"]]), ocProbesets())

+	AA <- A(callSet)

+	BB <- B(callSet)

+	for(j in 1:ncol(callSet)){

+		AA[snp.index, j] <- snprmaRes[["A"]][, j]

+		BB[snp.index, j] <- snprmaRes[["B"]][, j]

+	}

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

+	cnrmaRes <- cnrma2(A=AA,

 			   filenames=filenames,

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

 			   cdfName=cdfName,

@@ -261,32 +264,169 @@ genotype2 <- function(filenames,

 	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)

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

+	if(verbose) message("Leaving crlmmGT2")

 	open(tmp[["calls"]])

 	open(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)

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

+	GT <- snpCall(callSet)

+	GTP <- snpCallProbability(callSet)

+	for(j in 1:ncol(callSet)){

+		GT[snp.index, j] <- tmp[["calls"]][, j]

+		GTP[snp.index, j] <- tmp[["confs"]][, j]

+	}

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

 	callSet$gender <- tmp$gender

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

-	return(cnSet)

+	return(callSet)

+}

+

+

+

+genotype3 <- 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(!isPackageLoaded("ff")) stop("Must load package 'ff'")

+	if(!copynumber){

+		callSet <- crlmm2(filenames=filenames,

+				  cdfName=cdfName,

+				  mixtureSampleSize=mixtureSampleSize,

+				  eps=eps,

+				  verbose=verbose,

+				  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

+	##lM(callSet) <- initializeParamObject(list(featureNames(callSet), unique(protocolData(callSet)$batch)))

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

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

+##	snprmaRes <- snprma2(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")

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

+##	##batches <- splitIndicesByLength(1:nrow(snprmaRes[["A"]]), ocProbesets())

+##	for(j in 1:ncol(callSet)){

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

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

+##	}

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

+##	if(verbose) message("Entering crlmmGT2...")

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

+##	if(verbose) message("Leaving crlmmGT2")

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

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

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

+##	for(j in 1:ncol(callSet)){

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

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

+##	}

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

+##	callSet$gender <- tmp$gender

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

+##	return(callSet)

+	return(cnrmaRes)

 }

+

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

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

 ## For Illumina

@@ -758,6 +898,7 @@ predictGender <- function(res, cdfName="genomewidesnp6", SNRMin=5){

 }

+##replace with release/crlmm/R/cnrma-functions

 crlmmCopynumber <- function(object,

 			    chromosome=1:23,

 			    which.batches,

@@ -776,7 +917,8 @@ crlmmCopynumber <- function(object,

 			    MIN.PHI=2^3,

 			    THR.NU.PHI=TRUE,

 			    thresholdCopynumber=TRUE){

-	if(isPackageLoaded("ff")){

+	ffIsLoaded <- class(calls(object))[[1]] == "ff"

+	if(ffIsLoaded){

 		open(object)

 		open(object$SKW)

 		open(object$SNR)

@@ -787,7 +929,7 @@ crlmmCopynumber <- function(object,

 	stopifnot("isSnp" %in% fvarLabels(object))

 	##batch <- object$batch

 	batch <- batch(object)

-	if(isPackageLoaded("ff")) {

+	if(ffIsLoaded) {

 		open(object$SNR)

 		SNR <- object$SNR[, ]

 	} else SNR <-  object$SNR

@@ -844,13 +986,19 @@ crlmmCopynumber <- function(object,

 			fvarLabels(tmp)[13:14] <- paste(c("phiPrimeA", "phiPrimeB"), batchName, sep=".")

 			fvarLabels(tmp) <- gsub("_", ".", fvarLabels(tmp))

 			##fvarLabels(tmp) <- gsub("\\.[1-9]", "", fvarLabels(tmp))

-			fData(tmp) <- fData(tmp)[, fvarLabels(tmp) %in% names(physical(lM(object)))]

-			jj <- match(fvarLabels(tmp), names(lM(object)))

-			lM(object)[row.index, jj] <- fData(tmp)

-			##labels.asis <- fvarLabels(tmp)

-			##labels.asis <- gsub(paste("_", unique(tmp$batch), sep=""), paste(".", ii, sep=""), labels.asis)

-			##k <- match(labels.asis, colnames(lM(object)))

-			##lM(object)[row.index, k] <- fData(tmp)

+			if(ffIsLoaded){

+				fData(tmp) <- fData(tmp)[, fvarLabels(tmp) %in% names(physical(lM(object)))]

+				jj <- match(fvarLabels(tmp), names(lM(object)))

+				lM(object)[row.index, jj] <- fData(tmp)

+			} else {

+				nms <- paste(names(lM(object)), batchName, sep=".")

+				fData(tmp) <- fData(tmp)[, fvarLabels(tmp) %in% nms]

+				for(k in seq_along(fvarLabels(tmp))){

+					kk <- match(fvarLabels(tmp)[k], paste(names(lM(object)), batchName, sep="."))

+					column <- match(batchName, colnames(lM(object)[[k]]))

+					lM(object)[[k]][row.index, column] <- fData(tmp)[, k]

+				}

+			}			

 			rm(tmp); gc()

 			ii <- ii+1

 		}

@@ -879,12 +1027,15 @@ crlmmCopynumber2 <- function(object,

 	stopifnot("chromosome" %in% fvarLabels(object))

 	stopifnot("position" %in% fvarLabels(object))

 	stopifnot("isSnp" %in% fvarLabels(object))

+	lM(object) <- initializeParamObject(list(featureNames(object), unique(protocolData(object)$batch)))

 	batch <- batch(object)

-	XIndex.snps <- (1:nrow(object))[chromosome(object) == 23 & isSnp(object) & !is.na(chromosome(object))]

+	XIndex.snps <- which(chromosome(object) == 23 & isSnp(object) & !is.na(chromosome(object)))

 	##YIndex.snps <- (1:nrow(object))[chromosome(object) == 24 & isSnp(object)]

-	XIndex.nps <- (1:nrow(object))[chromosome(object) == 23 & !isSnp(object) & !is.na(chromosome(object))]

-	autosomeIndex.snps <- (1:nrow(object))[chromosome(object) < 23 & isSnp(object) & !is.na(chromosome(object))]

-	autosomeIndex.nps <- (1:nrow(object))[chromosome(object) < 23 & !isSnp(object) & !is.na(chromosome(object))]	

+	XIndex.nps <- which(chromosome(object) == 23 & !isSnp(object) & !is.na(chromosome(object)))

+	##autosomeIndex.snps <- (1:nrow(object))[chromosome(object) < 23 & isSnp(object) & !is.na(chromosome(object))]

+	autosomeIndex.snps <- which(chromosome(object) < 23 & isSnp(object) & !is.na(chromosome(object)))

+	autosomeIndex.nps <- which(chromosome(object) < 23 & !isSnp(object) & !is.na(chromosome(object)))

+	##autosomeIndex.nps <- (1:nrow(object))[chromosome(object) < 23 & !isSnp(object) & !is.na(chromosome(object))]	

 	## Do chromosome X in batches

 	Ns <- initializeBigMatrix("Ns", nrow(object), 5)

@@ -1017,8 +1168,8 @@ fit.lm1 <- function(idxBatch,

 	open(snpflags)

 	open(normal)

-	corr <- corrA.BB <- corrB.AA <- sig2B <- sig2A <- tau2B <- tau2A <- matrix(NA, length(snps), length(unique(batch(object))))

-	flags <- nuA <- nuB <- phiA <- phiB <- corr

+	corrAB <- corrBB <- corrAA <- sig2B <- sig2A <- tau2B <- tau2A <- matrix(NA, length(snps), length(unique(batch(object))))

+	flags <- nuA <- nuB <- phiA <- phiB <- corrAB

 	normal.snps <- normal[snps, ]

 	cB <- cA <- matrix(NA, length(snps), ncol(object))

@@ -1053,9 +1204,9 @@ fit.lm1 <- function(idxBatch,

 		tau2B[, J] <- shrink(taus[, 3, drop=FALSE], Ns[, 1], DF.PRIOR)

 		sig2B[, J] <- shrink(taus[, 1, drop=FALSE], Ns[, 3], DF.PRIOR)

-		corr[, J] <- corByGenotype(A=A, B=B, G=G, Ns=Ns, which.cluster=2, DF.PRIOR)

-		corrB.AA[, J] <- corByGenotype(A=A, B=B, G=G, Ns=Ns, which.cluster=1, DF.PRIOR)

-		corrA.BB[, J] <- corByGenotype(A=A, B=B, G=G, Ns=Ns, which.cluster=3, DF.PRIOR)

+		corrAB[, J] <- corByGenotype(A=A, B=B, G=G, Ns=Ns, which.cluster=2, DF.PRIOR)

+		corrAA[, J] <- corByGenotype(A=A, B=B, G=G, Ns=Ns, which.cluster=1, DF.PRIOR)

+		corrBB[, J] <- corByGenotype(A=A, B=B, G=G, Ns=Ns, which.cluster=3, DF.PRIOR)

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

 		## Impute sufficient statistics for unobserved genotypes (plate-specific)

@@ -1067,7 +1218,8 @@ fit.lm1 <- function(idxBatch,

 		size <- min(5000, length(index.complete))

 		if(size == 5000) index.complete <- sample(index.complete, 5000, replace=TRUE)

 		if(length(index.complete) < 200){

-			stop("fewer than 200 snps pass criteria for predicting the sufficient statistics")

+			warning("fewer than 200 snps pass criteria for predicting the sufficient statistics")

+			return()

 		}

 		index <- vector("list", 3)

 		index[[1]] <- which(Ns[, 1] == 0 & (Ns[, 2] >= MIN.OBS & Ns[, 3] >= MIN.OBS))

@@ -1105,6 +1257,7 @@ fit.lm1 <- function(idxBatch,

 		flags[, J] <- rowSums(Ns == 0) > 0

 		##flags[, J] <- index[[1]] | index[[2]] | index[[3]] | rowSums(

+		## we're regressing on the medians using the standard errors (hence the division by N) as weights

 		##formerly coefs()

 		Np <- Ns

 		Np[Np < 1] <- 1

@@ -1176,6 +1329,16 @@ fit.lm1 <- function(idxBatch,

 	tmp <- physical(lM(object))$phiB

 	tmp[snps, ] <- phiB

 	lM(object)$phiB <- tmp

+	tmp <- physical(lM(object))$corrAB

+	tmp[snps, ] <- corrAB

+	lM(object)$corrAB <- tmp

+	tmp <- physical(lM(object))$corrAA

+	tmp[snps, ] <- corrAA

+	lM(object)$corrAA <- tmp

+	tmp <- physical(lM(object))$corrBB

+	tmp[snps, ] <- corrBB

+	lM(object)$corrAB <- tmp

+	

 	lapply(assayData(object), close)

 	lapply(lM(object), close)

 	TRUE

@@ -1380,13 +1543,14 @@ fit.lm3 <- function(idxBatch,

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

 		## Impute sufficient statistics for unobserved genotypes (plate-specific)

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

-		index <- apply(Ns.F, 2, function(x, MIN.OBS) which(x > MIN.OBS), MIN.OBS)

+		index <- apply(Ns.F, 2, function(x, MIN.OBS) which(x >= MIN.OBS), MIN.OBS)

 		correct.orderA <- muA.F[, 1] > muA.F[, 3]

 		correct.orderB <- muB.F[, 3] > muB.F[, 1]

 		index.complete <- intersect(which(correct.orderA & correct.orderB), intersect(index[[1]], intersect(index[[2]], index[[3]])))

 		size <- min(5000, length(index.complete))

 		if(length(index.complete) < 200){

-			stop("fewer than 200 snps pass criteria for predicting the sufficient statistics")

+			warning("fewer than 200 snps pass criteria for predicting the sufficient statistics")

+			return()

 		}

 		if(size==5000) index.complete <- sample(index.complete, size)

 		index <- vector("list", 3)

@@ -1520,6 +1684,16 @@ fit.lm3 <- function(idxBatch,

 	tmp <- physical(lM(object))$phiPrimeB

 	tmp[snps, ] <- phiB2

 	lM(object)$phiPrimeB <- tmp

+	

+	tmp <- physical(lM(object))$corrAB

+	tmp[snps, ] <- corrAB

+	lM(object)$corrAB <- tmp

+	tmp <- physical(lM(object))$corrAA

+	tmp[snps, ] <- corrAA

+	lM(object)$corrAA <- tmp

+	tmp <- physical(lM(object))$corrBB

+	tmp[snps, ] <- corrBB

+	lM(object)$corrAB <- tmp	

 ##	lM(object)$tau2A[snps, ] <- tau2A

 ##	lM(object)$tau2B[snps, ] <- tau2B

 ##	lM(object)$sig2A[snps, ] <- sig2A

@@ -54,7 +54,6 @@ readIdatFiles <- function(sampleSheet=NULL,

                }

                pd = new("AnnotatedDataFrame", data = data.frame(Sample_ID=arrayNames))

        }

-       

        narrays = length(arrayNames)

        grnfiles = paste(arrayNames, fileExt$green, sep=sep)

        redfiles = paste(arrayNames, fileExt$red, sep=sep)

@@ -360,7 +359,6 @@ readIdatFiles2 <- function(sampleSheet=NULL,

 ## the readIDAT() and readBPM() functions below were provided by Keith Baggerly, 27/8/2008

 readIDAT <- function(idatFile){

-

   fileSize <- file.info(idatFile)$size

   tempCon <- file(idatFile,"rb")

@@ -91,18 +91,34 @@ setMethod("computeCopynumber", "CNSet",

 	}

 	object

 })

-

-

 setMethod("copyNumber", "CNSet", function(object){

 	I <- isSnp(object)

+	ffIsLoaded <- class(calls(object))[[1]]=="ff"

 	CA <- CA(object)

 	CB <- CB(object)

+	if(ffIsLoaded){

+		open(CA)

+		open(CB)

+		CA <- CA[,]

+		CB <- CB[,]

+	}

 	CN <- CA + CB

 	##For nonpolymorphic probes, CA is the total copy number

 	CN[!I, ] <- CA(object)[!I, ]

+	CN <- CN/100

 	CN

 })

+##setMethod("copyNumber", "CNSet", function(object){

+##	I <- isSnp(object)

+##	CA <- CA(object)

+##	CB <- CB(object)

+##	CN <- CA + CB

+##	##For nonpolymorphic probes, CA is the total copy number

+##	CN[!I, ] <- CA(object)[!I, ]

+##	CN

+##})

+

 setMethod("ellipse", "CNSet", function(x, copynumber, batch, ...){

 	ellipse.CNSet(x, copynumber, batch, ...)

new file mode 100644

Binary files /dev/null and b/data/sample.CNSetLM.rda differ

@@ -13,6 +13,7 @@

 \newcommand{\Rclass}[1]{{\textit{#1}}}

 \newcommand{\oligo}{\Rpackage{oligo }}

 \newcommand{\R}{\textsf{R}}

+\newcommand{\crlmm}{\Rpackage{crlmm}}

 \usepackage[margin=1in]{geometry}

 \begin{document}

@@ -22,8 +23,7 @@

 \maketitle

 <<setup, echo=FALSE, results=hide>>=

-options(continue=" ")

-options(prompt="R> ")

+options(continue=" ", width=70)

 @ 

 %\section{Estimating copy number}

@@ -43,168 +43,304 @@ options(prompt="R> ")

 CRLMM supports the following platforms:

 <<cdfname>>=

+library(oligoClasses)

 library(crlmm)

 crlmm:::validCdfNames()

+cdfName <- "genomewidesnp6"

 @ 

 \paragraph{Preprocess and genotype.}

-Provide the complete path for the filenames:

+In the following code chunk, we provide the complete path to the

+Affymetrix CEL files, assign a path to store the normalized

+intensities and genotype data, and define a 'batch' variable. The

+batch variable will later be useful when we estimate copy number.  We

+typically use the 96 well chemistry plate or the scan date of the

+array as a surrogate for batch. Adjusting by date or chemistry plate

+can be helpful for limiting the influence of batch effects.  

 <<celfiles>>=

 celFiles <- list.celfiles("/thumper/ctsa/snpmicroarray/hapmap/raw/affy/1m", full.names=TRUE, pattern=".CEL")

-## CEPH and Yoruban files

+outdir <- "/thumper/ctsa/snpmicroarray/rs/data/hapmap/1m/crlmm"

+if(!file.exists(outdir)) dir.create(outdir)

+## CEPH and Yoruban

 batch <- substr(basename(celFiles), 13, 13)

 celFiles <- celFiles[batch %in% c("C", "Y")]

 batch <- batch[batch %in% c("C", "Y")]

+stopifnot(length(batch) == length(celFiles))

 @ 

-While genotyping with crlmm can be performed using a small number of

-samples, copy number estimation requires at least 10 samples --

-preferably all of the samples that were processed together on the same

-plate.  Crlmm does not use a reference dataset when estimating model

-parameters because of large batch effects \citep{Scharpf2009}. The

-quantile normalization performed as part of the preprocessing of the raw

-data is insufficient for removing batch effects.  Processing a reference

-dataset, such as HapMap samples, along with the experimental data will

-not improve copy number estimation for the experimental dataset, and

-should not be used as a means to increase the sample size. Furthermore,

-processing a reference dataset without acknowledging that these samples

-were derived from a different batch can result in incorrect copy number

-estimates in both the experimental and reference datasets.  The

-appropriate way to acknowledge batch is to supply the batch name for

-each sample to be processed in the argument to \Robject{cnOptions}:

-

-<<celfiles>>=

-cnOpts <- cnOptions(cdfName="genomewidesnp6",

-		    outdir="/thumper/ctsa/beaty/scharpf/crlmmOut/hapmap",

-		    batch=batch)

-str(cnOpts)

-stopifnot(length(cnOpts$batch) == length(celFiles))

-names(cnOpts$batch) <- basename(celFiles)

+The preprocessing steps for copy number estimation includes quantile

+normalization of the polymorphic and nonpolymorphic markers and a

+summarization step whereby the quantile normalized intenities are

+summarized for each locus.  As the markers at polymorphic loci on the

+Affymetrix 6.0 platform are identical, we summarize the intensities by

+the median. For the nonpolymorphic markers, no summarization step is

+performs.  Next, the \crlmm{} package provides a genotype call and a

+confidence score at each polymorphic locus.  Unless the dataset is small

+(e.g., fewer than 50 samples), we suggest installing and loading the

+\R{} package \Rpackage{ff}.  The function \Rfunction{genotype} checks to

+see whether the \Rpackage{ff} is loaded.  If loaded, the normalized

+intensities and genotype are stored as \Robject{ff} objects on disk.

+Otherwise, the genotypes and normalized intensities are stored in

+matrices.  A word of caution: the \Rfunction{genotype} function requires

+a potentially large amount of RAM even when the \R{} package

+\Rpackage{ff} is loaded.  Users with large datasets may want to skip

+this part.  For the 180 HapMap samples processed in this vignette, the

+\Rfunction{genotype} function required 25MB RAM. The next two code

+chunks should be submitted as batch jobs.

+

+<<genotype>>=

+if(!exists("cnSet")){

+	if(!file.exists(file.path(outdir, "cnSet.rda"))){

+		##ops <- options(error=recover)

+		cnSet <- genotype(celFiles, cdfName=cdfName, copynumber=TRUE, batch=batch)

+		save(cnSet, file=file.path(outdir, "cnSet.rda"))

+	} else load(file.path(outdir, "cnSet.rda"))

+}

 @ 

-\noindent The next code chunk quantile normalizes the samples to a

-target reference distribution, uses the crlmm algorithm to genotype, and

-then estimates the copy number for each batch. Currently processing the

-180 HapMap cel files will require less than 20G of RAM. We are working

-on methods to reduce the memory footprint.

+A more memory efficient approach to preprocessing and genotyping is

+implemented in the \R{} function \Rfunction{genotype2}.  The arguments

+to this function are similar to \Rfunction{genotype} but requires

+explicit loading of the \R{} package \Rpackage{ff} prior to calling

+the function. The functions \Rfunction{ocProbesets} and

+\Rfunction{ocSamples} can be used to manage how many probesets and

+samples are to be loaded into memory and processed.  Using the default

+settings, the following code required 1.9Mb RAM to process 180 CEL

+files.

+

+<<genotype2>>=

+library(ff)

+ld.path <- file.path(outdir, "ffobjs")

+if(!file.exists(ld.path)) dir.create(ld.path)

+ldPath(ld.path)

+if(!exists("cnSet2")){

+	if(!file.exists(file.path(outdir, "cnSet2.rda"))){

+		cnSet2 <- genotype2(celFiles, cdfName=cdfName, copynumber=TRUE, batch=batch)

+		save(cnSet2, file=file.path(outdir, "cnSet2.rda"))

+	} else load(file.path(outdir, "cnSet2.rda"))

+}

+@ 

-<<preprocessAndGenotype, echo=TRUE>>=

-if(FALSE) crlmmCopynumber(celFiles, cnOpts) 

-load(file.path(cnOpts[["outdir"]], "cnSet_21.rda"))

+The objects returned by the \Rfunction{genotype} and

+\Rfunction{genotype2} differ in size as the former returns an object

+with matrices in the assay data slot, whereas the latter return an

+object with pointers to files on disk.  The functions \Rfunction{open}

+and \Rfunction{close} can be used to open and close the connections

+stored in the assay data of the \Robject{cnSet2} object, respectively.

+Subsetting the \Robject{ff} object pulls the data from disk into

+memory.  As subsetting operations pull the data from disk to memory,

+these operations must be performed with care.  In particular,

+subsetting the \Robject{cnSet2} will subset each assay data element

+and this can be slow.  The preferred approach would be to extract the

+assay data element that is needed, and then subset this function

+object as illustrated below for the genotype calls.

+

+<<check>>=

+class(calls(cnSet))

+dims(cnSet)

+class(calls(cnSet2))

+dims(cnSet2)

+print(object.size(cnSet), units="Mb")

+print(object.size(cnSet2), units="Mb")

+if(FALSE){

+	replicate(5, system.time(gt1 <- calls(cnSet)[, 1:50])[[1]])

+	open(calls(cnSet2))

+	replicate(5, system.time(gt2 <- calls(cnSet2)[, 1:50])[[1]])

+}

+gt1 <- calls(cnSet)[, 1:50]

+gt2 <- calls(cnSet2)[, 1:50]

+all.equal(gt1, gt2)

 @ 

-The following R objects are created from crlmmCopynumber:

+\noindent With \Robject{ff} objects the additional I/O time required

+for extracting data is substantial and will increase for larger

+datasets.  Patience is required.

+

+Note that for the Affymetrix 6.0 platform the assay data elements each

+have a row dimension corresponding to the total number of polymorphic

+and nonpolymorphic markers interrogated by the Affymetrix 6.0 platform.

+A consequence of keeping the rows of the assay data elements the same

+for all of the statistical summaries is that the matrix used to store

+genotype calls is larger than necessary.  Also, note the additional

+overhead of some operations when using \Robject{ff} objects.  For

+instance, the posterior probabilities for the CRLMM genotype calls are

+represented as integers. The accessor \Robject{snpCallProbability} can

+be used to access these confidence scores.  When stored as matrices,

+converting the integer representation back to the probability scale is

+straightforward as shown below.  However, for the \Robject{ff} objects

+we must first bring the data into memory. To bring all the scores into

+memory, one could use the function \Rfunction{[,]} but this could be

+slow and require a lot of RAM depending on the size of the dataset.

+Instead, we suggest pulling only the needed rows and columns from

+memory. In the following example, we convert the integer scores to

+probabilities for the CEPH samples.  As genotype confidence scores are

+not applicable to the nonpolymorphic markers, we extract only the

+polymorphic markers using the \Rfunction{isSnp} function.

+

+<<assayData>>=

+rows <- which(isSnp(cnSet))

+cols <- which(batch == "C")

+posterior.prob <- tryCatch(integerScoreToProbability(snpCallProbability),

+			    error=function(e) print("This will not work for an ff object."))

+posterior.prob1 <- integerScoreToProbability(snpCallProbability(cnSet)[rows, cols])

+posterior.prob2 <- integerScoreToProbability(snpCallProbability(cnSet2)[rows, cols])

+all.equal(posterior.prob1, posterior.prob2)

+@ 

-<<crlmmSetListObjects>>=

-fns <- list.files(cnOpts[["outdir"]], pattern="cnSet", full.name=TRUE)

-basename(fns)[1:5]

+Next, we obtain locus-level estimates of copy number by fitting a linear

+model to each SNP. A variable named 'batch' must be indicated in the

+\Robject{protocolData} of the \Robject{cnSet} object. As the inverse

+variance of the within-genotype normalized intensities are used as

+weights in the linear model (and hence the design matrix in the

+regression changes for each SNP), the time is linear with the number of

+markers on the array. Copy number estimation at nonpolymorphic markers

+and polymorphic markers with unobserved genotypes is more difficult. We

+refer the interested reader to the technical report \citep{Scharpf2009}.

+Again, we peform the copy number estimation using the matrix version and

+the ff version in parallel and encourage users with large datasets to

+explore the latter. As with the preprocessing and genotyping, the

+following code should be submitted as part of the batch job as it is too

+slow for interactive analysis.

+

+<<cn, echo=FALSE, eval=FALSE>>=

+##Matrix format

+cnSet <- crlmmCopynumber(cnSet)

 @ 

+<<save.cnset, echo=FALSE, eval=FALSE>>=

+save(cnSet, file=file.path(outdir, "cnSet.rda"))

+@ 

-The above algorithm for estimating copy number is predicated on the

-assumption that most samples within a batch have copy number 2 at any

-given locus.  For common copy number variants, this assumption may not

-hold.  An additional iteration using a bias correction provides

-additional robustness to this assumption.  Set the \Robject{bias.adj}

-argument to \Robject{TRUE}:

+<<cn.ff, echo=FALSE, eval=FALSE>>=

+rm(cnSet); gc()

+ocProbesets(75e3)

+##ff objects

+system.time(cnSet2 <- crlmmCopynumber2(cnSet2))

+save(cnSet2, file=file.path(outdir, "cnSet2.rda"))

+@ 

-<<biasAdjustment, eval=FALSE>>=

-cnOpts[["bias.adj"]] <- TRUE

-if(FALSE) crlmmCopynumber(celFiles, cnOpts)

+<<timings, eval=FALSE>>=

+tryCatch(print(paste("Time for matrix version:", time1)), error=function(e) print("timing for CN estimation not available"))