@@ -400,4 +400,8 @@ is stored in the assayData slot.

       - B.rda

       - cn_rmaResult.rda (SKW)

+2010-02-21 R. Scharpf - committed version 1.5.26

+ 

+** additional support for ff in crlmm and computeCopynumber functions

+

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

-Date: 2010-02-20

+Version: 1.5.26

+Date: 2010-02-21

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

 Maintainer: Benilton S Carvalho <bcarvalh@jhsph.edu>, 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

@@ -11,7 +11,8 @@ Depends: R (>= 2.11.0), methods, Biobase (>= 2.7.2), oligoClasses (>= 1.9.28)

 Imports: affyio (>= 1.15.2), preprocessCore, utils, stats, genefilter, splines, mvtnorm, ellipse, SNPchip

 Enhances: ff

 Suggests: hapmapsnp5, hapmapsnp6, genomewidesnp5Crlmm (>= 1.0.2),genomewidesnp6Crlmm (>= 1.0.2), snpMatrix, metaArray

-Collate: AllGenerics.R

+Collate: AllClasses.R

+	 AllGenerics.R

 	 methods-CNSet.R

 	 methods-eSet.R

          methods-SnpSuperSet.R

@@ -242,45 +242,47 @@ predictGender <- function(res, cdfName="genomewidesnp6", SNRMin=5){

-crlmmCopynumber <- function(filenames, cnOptions, object, ...){

-	cnOpts <- cnOptions

-	batchSize <- cnOptions[["crlmmBatchSize"]]

-	batch <- cnOptions[["batch"]]

-	outdir <- cnOptions[["outdir"]]

-	##if necessary, split plates in separate batches

-	if(length(filenames) > batchSize){

-		L <- length(table(batch))

-		message("Processing the files in")

-		batchList <- split(names(table(batch)), rep(1:L, each=10, length.out=L))

-		for(i in seq(along=batchList)){

-			index <- as.character(batch) %in% as.character(batchList[[i]])

-			fns <- filenames[index]

-			cnOpts[["batch"]] <- cnOptions[["batch"]][index]

-			cnOpts[["AFile"]] <- paste(outdir, "/", paste(cnOptions[["AFile"]], i, sep="_crlmmBatch"), ".rda", sep="")

-			cnOpts[["BFile"]] <- paste(outdir, "/", paste(cnOptions[["BFile"]], i, sep="_crlmmBatch"), ".rda", sep="")

-			cnOpts[["callsFile"]] <- paste(outdir, "/", paste(cnOptions[["callsFile"]], i, sep="_crlmmBatch"), ".rda", sep="")

-			cnOpts[["confsFile"]] <- paste(outdir, "/", paste(cnOptions[["confsFile"]], i, sep="_crlmmBatch"), ".rda", sep="")

-			cnOpts[["snprmaFile"]] <- paste(outdir, "/", paste(cnOptions[["snprmaFile"]], i, sep="_crlmmBatch"), ".rda", sep="")

-			cnOpts[["cnrmaFile"]] <- paste(outdir, "/", paste(cnOptions[["cnrmaFile"]], i, sep="_crlmmBatch"), ".rda", sep="")

-			cnOpts[["cnFile"]] <- file.path(outdir, paste(cnOptions[["cnFile"]], i, sep="_crlmmBatch"))

-			cnOpts[["rgFile"]] <- file.path(outdir, paste(cnOptions[["rgFile"]], i, sep="_crlmmBatch"))

-			crlmmCopynumber.crlmmBatch(fns, cnOpts, object, ...)

-		}

-	} else {

-		cnOpts[["AFile"]] <- paste(outdir, "/", cnOptions[["AFile"]], ".rda", sep="")

-		cnOpts[["BFile"]] <- paste(outdir, "/", cnOptions[["BFile"]], ".rda", sep="")

-		cnOpts[["callsFile"]] <- paste(outdir, "/", cnOptions[["callsFile"]], ".rda", sep="")

-		cnOpts[["confsFile"]] <- paste(outdir, "/", cnOptions[["confsFile"]], ".rda", sep="")

-		cnOpts[["snprmaFile"]] <- paste(outdir, "/",cnOptions[["snprmaFile"]], ".rda", sep="")

-		cnOpts[["cnrmaFile"]] <- paste(outdir, "/", cnOptions[["cnrmaFile"]], ".rda", sep="")

-		cnOpts[["cnFile"]] <- file.path(outdir, cnOptions[["cnFile"]])

-		cnOpts[["rgFile"]] <- file.path(outdir, cnOptions[["rgFile"]])

-		crlmmCopynumber.crlmmBatch(filenames, cnOpts, object, ...)

-	}

-}

+##crlmmCopynumber <- function(filenames, cnOptions, object, ...){

+##	cnOpts <- cnOptions

+##	batchSize <- cnOptions[["crlmmBatchSize"]]

+##	batch <- cnOptions[["batch"]]

+##	outdir <- cnOptions[["outdir"]]

+##	##if necessary, split plates in separate batches

+##	if(length(filenames) > batchSize){

+##		L <- length(table(batch))

+##		message("Processing the files in")

+##		batchList <- split(names(table(batch)), rep(1:L, each=10, length.out=L))

+##		for(i in seq(along=batchList)){

+##			index <- as.character(batch) %in% as.character(batchList[[i]])

+##			fns <- filenames[index]

+##			cnOpts[["batch"]] <- cnOptions[["batch"]][index]

+##			cnOpts[["AFile"]] <- paste(outdir, "/", paste(cnOptions[["AFile"]], i, sep="_crlmmBatch"), ".rda", sep="")

+##			cnOpts[["BFile"]] <- paste(outdir, "/", paste(cnOptions[["BFile"]], i, sep="_crlmmBatch"), ".rda", sep="")

+##			cnOpts[["CAFile"]] <- paste(outdir, "/", paste(cnOptions[["CAFile"]], i, sep="_crlmmBatch"), ".rda", sep="")

+##g			cnOpts[["CBFile"]] <- paste(outdir, "/", paste(cnOptions[["CBFile"]], i, sep="_crlmmBatch"), ".rda", sep="")			

+##			cnOpts[["callsFile"]] <- paste(outdir, "/", paste(cnOptions[["callsFile"]], i, sep="_crlmmBatch"), ".rda", sep="")

+##			cnOpts[["confsFile"]] <- paste(outdir, "/", paste(cnOptions[["confsFile"]], i, sep="_crlmmBatch"), ".rda", sep="")

+##			cnOpts[["snprmaFile"]] <- paste(outdir, "/", paste(cnOptions[["snprmaFile"]], i, sep="_crlmmBatch"), ".rda", sep="")

+##			cnOpts[["cnrmaFile"]] <- paste(outdir, "/", paste(cnOptions[["cnrmaFile"]], i, sep="_crlmmBatch"), ".rda", sep="")

+##			cnOpts[["cnFile"]] <- file.path(outdir, paste(cnOptions[["cnFile"]], i, sep="_crlmmBatch"))

+##			cnOpts[["rgFile"]] <- file.path(outdir, paste(cnOptions[["rgFile"]], i, sep="_crlmmBatch"))

+##			crlmmCopynumber.crlmmBatch(fns, cnOpts, object, ...)

+##		}

+##	} else {

+##		cnOpts[["AFile"]] <- paste(outdir, "/", cnOptions[["AFile"]], ".rda", sep="")

+##		cnOpts[["BFile"]] <- paste(outdir, "/", cnOptions[["BFile"]], ".rda", sep="")

+##		cnOpts[["callsFile"]] <- paste(outdir, "/", cnOptions[["callsFile"]], ".rda", sep="")

+##		cnOpts[["confsFile"]] <- paste(outdir, "/", cnOptions[["confsFile"]], ".rda", sep="")

+##		cnOpts[["snprmaFile"]] <- paste(outdir, "/",cnOptions[["snprmaFile"]], ".rda", sep="")

+##		cnOpts[["cnrmaFile"]] <- paste(outdir, "/", cnOptions[["cnrmaFile"]], ".rda", sep="")

+##		cnOpts[["cnFile"]] <- file.path(outdir, cnOptions[["cnFile"]])

+##		cnOpts[["rgFile"]] <- file.path(outdir, cnOptions[["rgFile"]])

+##		crlmmCopynumber.crlmmBatch(filenames, cnOpts, object, ...)

+##	}

+##}

-crlmmCopynumber.crlmmBatch <- function(filenames, cnOptions, object, ...){

+crlmmCopynumber <- function(filenames, cnOptions, object, ...){

 	set.seed(cnOptions[["seed"]])  ##for reproducibility

 	if(!missing(object)){

 		stopifnot(class(object) == "CNSet")

@@ -292,24 +294,30 @@ crlmmCopynumber.crlmmBatch <- function(filenames, cnOptions, object, ...){

 		##crlmmResults <- crlmmWrapper(filenames, cnOptions, ...)

 		crlmmWrapper(filenames, cnOptions, ...)

 	}

+	verbose <- cnOptions[["verbose"]]

 	load.it <- cnOptions[["load.it"]]

 	outdir <- cnOptions[["outdir"]]

 	snprmaFile <- cnOptions[["snprmaFile"]]

 	cnFile <- cnOptions[["cnFile"]]

 	cnrmaFile <- cnOptions[["cnrmaFile"]]

+	snprmaFile <- cnOptions[["snprmaFile"]]

 	callsFile <- cnOptions[["callsFile"]]

 	confsFile <- cnOptions[["confsFile"]]

 	AFile <- cnOptions[["AFile"]]

 	BFile <- cnOptions[["BFile"]]

+	CAFile <- cnOptions[["CAFile"]]

+	CBFile <- cnOptions[["CBFile"]]

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

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

 	snpProbes <- get("snpProbes")

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

 	cnProbes <- get("cnProbes")	

 	k <- grep("chr", colnames(snpProbes))

-	message("Attaching quantile-normalized intensities...")

+	if(verbose) message("Loading quantile-normalized intensities...")

 	##cwd <- getwd()

 	##setwd(dirname(snprmaFile[1]))

+	load(snprmaFile)

+	res <- get("res")

 	load(AFile)

 	if(isPackageLoaded("ff")) open(A)

 	load(BFile)

@@ -321,28 +329,27 @@ crlmmCopynumber.crlmmBatch <- function(filenames, cnOptions, object, ...){

 	load(confsFile)

 	confs <- get("confs")

 	if(isPackageLoaded("ff")) open(confs)	

-	fns <- rownames(A)

+	##fns <- rownames(A)

+	fns <- c(res[["gns"]], res[["cnnames"]])

+	if(length(fns) != nrow(A)) stop("check featurenames. fns should be the rownames of A")

 	nr <- nrow(A); nc <- ncol(A)

 	if(isPackageLoaded("ff")){

-		##initialize CA, CB

-		if(!file.exists(file.path(outdir, "CA.rda"))){

-			CA <- ff(dim=c(nr, nc), vmode="integer", finalizer="close", dimnames=dimnames(A))

-

-		} else{

-			load(file.path(outdir, "CA.rda"))

-			if(!all(dim(CA) == dim(A)))

-				CA <- ff(dim=c(nr, nc), vmode="integer", finalizer="close", dimnames=dimnames(A), overwrite=TRUE)				

-			open(CA)

-		}

-		if(!file.exists(file.path(outdir, "CB.rda"))){		

-			CB <- ff(dim=c(nr, nc), vmode="integer", finalizer="close", dimnames=dimnames(A))

-		} else{

-			load(file.path(outdir, "CB.rda"))

-			if(!all(dim(CB) == dim(A)))

-				CB <- ff(dim=c(nr, nc), vmode="integer", finalizer="close", dimnames=dimnames(A), overwrite=TRUE)							

-			open(CB)

-		}

-	}

+		if(file.exists(CAFile)){

+			load(CAFile)

+			if(is.null(dim(CA)) | !all(dim(CA) == dim(A))) {

+				unlink(CAFile)

+				CA <- initializeBigMatrix(nr, nc)

+			}

+		} else CA <- initializeBigMatrix(nr, nc)

+		if(file.exists(CBFile)){

+			load(CBFile)

+			if(is.null(dim(CB)) | !all(dim(CB) == dim(A))) {			     

+				CB <- initializeBigMatrix(nr, nc)

+			}

+		} else CB <- initializeBigMatrix(nr, nc)

+		open(CA)

+		open(CB)

+	} 

 	##cnFile <- cnOptions[["cnFile"]]

 	chromosome <- cnOptions[["chromosome"]]

 	SNRmin <- cnOptions[["SNRmin"]]

@@ -353,29 +360,47 @@ crlmmCopynumber.crlmmBatch <- function(filenames, cnOptions, object, ...){

 			    data=scanDates,

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

 			    row.names=colnames(scanDates)))

-	load(file.path(cnOptions[["outdir"]], "sampleStats.rda"))

-	sampleStats <- data.frame(sampleStats, batch=cnOptions[["batch"]])

+	##load(file.path(cnOptions[["outdir"]], "sampleStats.rda"))

+	sampleStats <- data.frame(SKW=res$SKW,

+				  SNR=res$SNR,

+				  mixtureParams=res$mixtureParams,

+				  gender=res$gender,

+				  batch=cnOptions[["batch"]])

+	##sampleStats, batch=cnOptions[["batch"]])

 	pD <- new("AnnotatedDataFrame",

 		  data=sampleStats,

 		  varMetadata=data.frame(labelDescription=colnames(sampleStats)))

 	batch <- cnOptions[["batch"]]

+	if(isPackageLoaded("ff")){

+##		featureDataFF <- ffdf(position=ff(rep(0L, length(fns)), vmode="integer", finalizer="close"),

+##				      chromosome=ff(rep(0L, length(fns)), vmode="integer", finalizer="close"),

+##				      isSnp=ff(rep(0L, length(fns)), vmode="integer", finalizer="close"),

+##				      row.names=fns)

+##		fD <- vector("list", length(chromosome))

+		featureDataFF <- ff(dim=c(nrow(A), 17*length(unique(batch))+3), vmode="double", finalizer="close",

+				    overwrite=TRUE)

+	}

 	for(CHR in chromosome){

 		snps <- rownames(snpProbes)[snpProbes[, k] == CHR]

 		cns <- rownames(cnProbes)[cnProbes[, k] == CHR]		

-		index.snps <- match(snps, rownames(A))

-		index.cn <- match(cns, rownames(A))

+		index.snps <- match(snps, fns)

+		index.cn <- match(cns, fns)

 		row.index <- c(index.snps, index.cn)

 		cat("Chromosome ", CHR, "\n")

+		if(isPackageLoaded("ff")){		

+			fD.batch <- vector("list", length(unique(batch)))

+		}

 		for(i in seq(along=unique(batch))){

 			PLATE <- unique(batch)[i]

+			message("Plate: ", PLATE)

 			sample.index <- which(batch==PLATE)

 			##cnOptions[["batch"]] <- cnOptions[["batch"]][snpI[["SNR"]]  >= SNRmin]

 			if(isPackageLoaded("ff")){

 				ca <- CA[row.index, sample.index]

 				cb <- CB[row.index, sample.index]

 			} else{

-				dns <- dimnames(A[row.index, sample.index])

-				cb <- ca <- matrix(NA, nr=length(row.index), nc=length(sample.index), dimnames=dns)

+				##dns <- dimnames(A[row.index, sample.index])

+				cb <- ca <- matrix(NA, nr=length(row.index), nc=length(sample.index))##, dimnames=dns)

 			}

 			cnSet <- new("CNSet",

 				     alleleA=A[row.index, sample.index],

@@ -386,52 +411,64 @@ crlmmCopynumber.crlmmBatch <- function(filenames, cnOptions, object, ...){

 				     CB=cb,

 				     featureData=annotatedDataFrameFrom(A[row.index, sample.index], byrow=TRUE),

 				     phenoData=pD[sample.index, ],

-				     protocolData=protocoldata[sample.index, ],

-				     annotation=cnOptions[["cdfName"]])			

-			if(any(cnSet$SNR > SNRmin)){

-				message(paste("Excluding samples with SNR < ", SNRmin))

-				cnSet <- cnSet[, cnSet$SNR >= SNRmin]

-			}

-			##rm(snpI, NP, gtSet); gc()

+				     protocolData=protocoldata[sample.index, ])

+			##Verify this is correct

+			annotation(cnSet) <- cnOptions[["cdfName"]]

+			featureNames(cnSet) <- fns[row.index]

+			##add chromosome, position, isSnp

+			cnSet <- annotate(cnSet)

+##			if(any(cnSet$SNR > SNRmin)){

+##				if(CHR == chromosome[1]) message(paste("Excluding samples with SNR < ", SNRmin))

+##				cnSet <- cnSet[, cnSet$SNR >= SNRmin]

+##			}

 			featureData(cnSet) <- lm.parameters(cnSet, cnOptions)

-		} 

-		if(CHR != 24){

-			## 60G for 1239 files

-			cnSet <- computeCopynumber(cnSet, cnOptions)

-			if(isPackageLoaded("ff")){

-				CA[row.index, sample.index] <- cnSet@assayData[["CA"]]

-				CB[row.index, sample.index] <- cnSet@assayData[["CB"]]

-			}

-		} else{

-			message("Copy number estimates not available for chromosome Y.  Saving only the 'callSet' object for this chromosome")

-			alleleSet <- cnSet

-			rm(cnSet) ; gc()

-			save(alleleSet, file=file.path(cnOptions[["outdir"]], paste("alleleSet_", PLATE, "_", CHR, ".rda", sep="")))

-			rm(alleleSet); gc()

-			next()

-		}

-		if(length(chromosome) == 1){

-			if(cnOptions[["save.cnset"]]){

-				save(cnSet, file=paste(cnFile, "_", PLATE, "_", CHR, ".rda", sep=""))				

-			}

-		}

-		if(length(chromosome) > 1){

-			save(cnSet, file=paste(cnFile, "_", PLATE, "_", CHR, ".rda", sep=""))

-			##save(cnSet, file=file.path(cnOptions[["outdir"]], paste("cnSet_", CHR, ".rda", sep="")))

-			rm(cnSet); gc()

-		} else {

-			if(isPackageLoaded("ff")) {

-				close(CA); close(CB)

-				save(CA, file=file.path(cnOptions[["outdir"]], "CA.rda"))

-				save(CB, file=file.path(cnOptions[["outdir"]], "CB.rda"))		

-			}			

-			return(cnSet)

-		}

-	}

+			if(CHR < 24){

+				cnSet <- computeCopynumber(cnSet, cnOptions)

+				if(isPackageLoaded("ff")){

+					if(i == 1) {  ## first batch

+						featureDataFF[row.index, 1:20] <- as.matrix(fData(cnSet))

+						fcol <- 20

+					} else {

+						featureDataFF[row.index, (fcol+1):(fcol+17)]

+						fcol <- fcol+17

+						##remove redundant columns

+					}

+					CA[row.index, sample.index] <- cnSet@assayData[["CA"]]

+					CB[row.index, sample.index] <- cnSet@assayData[["CB"]]

+					next() ## next batch

+				} else {

+					save(cnSet, file=paste(cnFile, "_", PLATE, "_", CHR, ".rda", sep=""))

+				}

+			} else break()  ## break out of the batch loop.  Go to next chromosome

+			##		if(length(chromosome) == 1 & length(unique(batch)) == 1){

+			##			if(cnOptions[["save.cnset"]] & !isLoadedPackage("ff"))

+			##				save(cnSet, file=paste(cnFile, "_", PLATE, "_", CHR, ".rda", sep=""))

+			##			return(cnSet)

+			##		} else {

+			##			##either multiple chromosome or multiple batches...save to file if ff is not loaded

+			##			if(!isPackageLoaded("ff")){

+			##				save(cnSet, file=paste(cnFile, "_", PLATE, "_", CHR, ".rda", sep=""))

+			##			}

+			##			rm(cnSet); gc()

+			##	}

+		} ## end of batch loop

+	} ## end of chromosome loop

 	if(isPackageLoaded("ff")) {

+		close(A); close(B)

+		close(featureDataFF)

 		close(CA); close(CB)

-		save(CA, file=file.path(cnOptions[["outdir"]], "CA.rda"))

-		save(CB, file=file.path(cnOptions[["outdir"]], "CB.rda"))		

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

+		save(CA, file=CAFile)

+		save(CB, file=CBFile)

+		close(calls); close(confs)

+		ffSet <- new("FFSet",

+			     A=A,

+			     B=B,

+			     snpCall=calls,

+			     callProbability=confs,

+			     CA=CA,

+			     CB=CB)

+		return(ffSet)

 	}

 }

@@ -524,6 +561,7 @@ loadIlluminaCnrma <- function(){

 }

 crlmmWrapper <- function(filenames, cnOptions, ...){

+	crlmmBatchSize <- cnOptions[["crlmmBatchSize"]]

 	cdfName <- cnOptions[["cdfName"]]

 	load.it <- cnOptions[["load.it"]]

 	save.it <- cnOptions[["save.it"]]

@@ -548,25 +586,26 @@ crlmmWrapper <- function(filenames, cnOptions, ...){

 	}

 	if(platform == "illumina") RG <- loadIlluminaRG(rgFile, callsFile, load.it, save.it)

 	if(!(file.exists(dirname(callsFile)))) stop(dirname(callsFile), " does not exist.")

-	if(!(file.exists(dirname(snprmaFile[1])))) stop(dirname(snprmaFile[1]), " does not exist.")

-	if(load.it){

-		if(!file.exists(callsFile)){

-			message("load.it is TRUE, but ", callsFile, " does not exist.  Rerunning the genotype calling algorithm") 

-			##load.it <- FALSE

-		}

-	}

+	if(!(file.exists(dirname(snprmaFile)))) stop(dirname(snprmaFile), " does not exist.")

+##	if(load.it){

+##		if(!file.exists(callsFile)){

+##			message("load.it is TRUE, but ", callsFile, " does not exist.  Rerunning the genotype calling algorithm") 

+##			##load.it <- FALSE

+##		}

+##	}

 	if(platform == "affymetrix"){

-		if(!file.exists(callsFile) | !load.it){

+##		if(!file.exists(callsFile) | !load.it){

 			crlmm(filenames=filenames,

 			      cdfName=cdfName,

 			      save.it=TRUE,

 			      load.it=load.it,

-			      intensityFile=snprmaFile,

+			      snprmaFile=snprmaFile,

 			      callsFile=callsFile,

 			      confsFile=confsFile,

 			      AFile=AFile,

-			      BFile=BFile)

-		}

+			      BFile=BFile,

+			      crlmmBatchSize=crlmmBatchSize)

+##		}

 	}

 	gc()

 	if(platform == "illumina") {

@@ -756,18 +795,20 @@ thresholdCopynumber <- function(object){

 cnOptions <- function(outdir="./",

 		      cdfName,

-		      AFile="A",

-		      BFile="B",

-		      callsFile="genotypes",

-		      rgFile="rgFile",

+		      AFile="A.rda",

+		      BFile="B.rda",

+		      CAFile="CA.rda",

+		      CBFile="CB.rda",

+		      callsFile="genotypes.rda",

+		      rgFile="rgFile.rda",

 		      cnFile="cnSet",

-		      cnrmaFile="cn_rmaResult",

+		      cnrmaFile="cn_rmaResult.rda",

 		      ##npBinary="NP.bin",  ##name of file.backed.big.matrix

 		      ##npDesc="NP.desc",

-		      snprmaFile="snp_rmaResult",

+		      snprmaFile="snp_rmaResult.rda",

 		      ##callsFile="calls.desc",

 		      ##confsFile="confsFile.desc",

-		      confsFile="confsFile",

+		      confsFile="confs.rda",

 		      save.it=TRUE,

 		      save.cnset=TRUE,

 		      load.it=TRUE,

@@ -790,15 +831,16 @@ cnOptions <- function(outdir="./",

 		      thresholdCopynumber=TRUE,

 		      unlink=TRUE,

 		      use.poe=FALSE,

-		      crlmmBatchSize=1100, ## this is about right for Affy 6.0 platform 1.8e6 * 1100 < 2^31

+		      crlmmBatchSize=1000,

 		      ...){

-	if(isPackageLoaded("ff")){

-		AFile=paste(AFile, "ff", sep="_")

-		BFile=paste(BFile, "ff", sep="_")

-		confsFile=paste(confsFile, "ff", sep="_")

-		callsFile=paste(callsFile, "ff", sep="_")

-		cnFile=paste(cnFile, "ff", sep="")

-	}

+##	if(isPackageLoaded("ff")){

+##		AFile=paste(AFile, "ff", sep="_")

+##		BFile=paste(BFile, "ff", sep="_")

+##		confsFile=paste(confsFile, "ff", sep="_")

+##		callsFile=paste(callsFile, "ff", sep="_")

+##		cnFile=paste(cnFile, "ff", sep="")

+##	}

+	if(length(batch) > 200) warning("This job may require a lot of RAM. Consider using the ff package in conjuction with crlmm, as described in the copynumber_ff vignette...")

 	if(use.poe) require(metaArray)

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

 	if(!file.exists(outdir)){

@@ -816,14 +858,16 @@ cnOptions <- function(outdir="./",

 		stop("must specify batch -- should be the same length as the number of files")

 	list(outdir=outdir,

 	     cdfName=cdfName,

-	     callsFile=callsFile,

-	     confsFile=confsFile,

-	     AFile=AFile,

-	     BFile=BFile,

-	     snprmaFile=snprmaFile,

-	     cnrmaFile=cnrmaFile,

-	     rgFile=rgFile,

-	     cnFile=cnFile,

+	     callsFile=file.path(outdir, callsFile),

+	     confsFile=file.path(outdir, confsFile),

+	     AFile=file.path(outdir, AFile),

+	     BFile=file.path(outdir, BFile),

+	     CAFile=file.path(outdir, CAFile),

+	     CBFile=file.path(outdir, CBFile),

+	     snprmaFile=file.path(outdir, snprmaFile),

+	     cnrmaFile=file.path(outdir, cnrmaFile),

+	     rgFile=file.path(outdir, rgFile),

+	     cnFile=file.path(outdir, cnFile),

 	     save.it=save.it,

 	     save.cnset=save.cnset,

 	     load.it=load.it,

@@ -888,8 +932,11 @@ lm.parameters <- function(object, cnOptions){

 }

 nonpolymorphic <- function(object, cnOptions, tmp.objects){

+	chromosome <- cnOptions[["chromosome"]]

 	batch <- unique(object$batch)

 	CHR <- unique(chromosome(object))

+	verbose <- cnOptions[["verbose"]]

+	if(CHR != chromosome[1]) verbose <- FALSE

 	goodSnps <- function(object, PHI.THR, tmp.objects, nAA.THR, nBB.THR){

 		Ns <- tmp.objects[["Ns"]]

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

@@ -1015,10 +1062,9 @@ nonpolymorphic <- function(object, cnOptions, tmp.objects){

 		THR.NU.PHI <- cnOptions$THR.NU.PHI

 		if(THR.NU.PHI){

-			verbose <- cnOptions$verbose

 			##Assign values to object

 			object <- pr(object, "nuA", batch, nuA)

-			object <- pr(object, "phiA", batch, phiA)			

+			object <- pr(object, "phiA", batch, phiA)

 			if(verbose) message("Thresholding nu and phi")

 			object <- thresholdModelParams(object, cnOptions)

 		} else {

@@ -1036,7 +1082,6 @@ nonpolymorphic <- function(object, cnOptions, tmp.objects){

 		THR.NU.PHI <- cnOptions$THR.NU.PHI

 		if(THR.NU.PHI){

-			verbose <- cnOptions$verbose

 			##Assign values to object

 			object <- pr(object, "nuA", batch, nuA)

 			object <- pr(object, "phiA", batch, phiA)			

@@ -1738,6 +1783,7 @@ computeCopynumber.CNSet <- function(object, cnOptions){

 	if(length(batch) != ncol(object)) stop("Batch must be the same length as the number of samples")

 	MIN.SAMPLES <- cnOptions$MIN.SAMPLES

 	verbose <- cnOptions$verbose

+	if(CHR != cnOptions[["chromosome"]][1]) verbose <- FALSE

 	use.poe <- cnOptions[["use.poe"]]

 	for(i in seq(along=unique(batch))){

 		PLATE <- unique(batch)[i]

@@ -1754,11 +1800,11 @@ computeCopynumber.CNSet <- function(object, cnOptions){

 			cnOptions$bias.adj <- bias.adj <- FALSE

 		}

 		if(bias.adj){

-			if(verbose) message("Dropping samples with low posterior prob. of normal copy number (samples dropped is locus-specific)")

+			if(verbose & i == 1) message("Dropping samples with low posterior prob. of normal copy number (samples dropped is locus-specific)")

 			if(!use.poe)

 				tmp.objects <- biasAdjNP(object.batch, cnOptions, tmp.objects)

 			tmp.objects <- biasAdj(object.batch, cnOptions, tmp.objects)

-			if(verbose) message("Recomputing location and scale parameters")

+			if(verbose & i == 1) message("Recomputing location and scale parameters")

 		}

 		##update tmp.objects

 		tmp.objects <- withinGenotypeMoments(object.batch,

@@ -1773,14 +1819,13 @@ computeCopynumber.CNSet <- function(object, cnOptions){

 		##nuA=getParam(object.batch, "nuA", PLATE)

 		THR.NU.PHI <- cnOptions$THR.NU.PHI

 		if(THR.NU.PHI){

-			verbose <- cnOptions$verbose

-			if(verbose) message("Thresholding nu and phi")

+			if(verbose & i == 1) message("Thresholding nu and phi")

 			object.batch <- thresholdModelParams(object.batch, cnOptions)

 		}		

-		if(verbose) message("\nAllele specific copy number")	

+		if(verbose & i == 1) message("\nAllele specific copy number")	

 		object.batch <- polymorphic(object.batch, cnOptions, tmp.objects)

 		if(any(!isSnp(object))){ ## there are nonpolymorphic probes

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

+			if(verbose & i == 1) message("\nCopy number for nonpolymorphic probes...")

 			if(!use.poe){

 				object.batch <- nonpolymorphic(object.batch, cnOptions, tmp.objects)

 			} else {

@@ -1812,7 +1857,7 @@ computeCopynumber.CNSet <- function(object, cnOptions){

 		object <- pr(object, "corrB.AA", PLATE, getParam(object.batch, "corrB.AA", PLATE))		

 		rm(object.batch, tmp.objects); gc();

 	}

-	object <- object[order(chromosome(object), position(object)), ]

+	##object <- object[order(chromosome(object), position(object)), ]

 	if(cnOptions[["thresholdCopynumber"]]){

 		object <- thresholdCopynumber(object)

 	}

@@ -1,19 +1,23 @@

 crlmm <- function(filenames, row.names=TRUE, col.names=TRUE,

                   probs=c(1/3, 1/3, 1/3), DF=6, SNRMin=5, gender=NULL,

-                  save.it=FALSE, load.it=FALSE, intensityFile, callsFile, confsFile,

-		  AFile, BFile,

+                  save.it=FALSE, load.it=FALSE,

+		  snprmaFile="snp_rmaResult.rda",

+		  callsFile="genotypes.rda",

+		  confsFile="confs.rda",

+		  AFile="A.rda",

+		  BFile="B.rda",

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

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

-                  returnParams=FALSE, badSNP=0.7){

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

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

+                  returnParams=FALSE, badSNP=0.7,

+		  crlmmBatchSize=1000){

+	BS <- crlmmBatchSize

+	if(load.it & file.exists(snprmaFile) & file.exists(callsFile)) return() ##nothing to do

+	if (load.it & !file.exists(snprmaFile)){

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

+		message("'snprmaFile' does not exist and you chose to load.it.  Rerunning snprma...")

+		load.it <- FALSE

+	} 

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

-	if (!missing(intensityFile))

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

-			load.it <- FALSE

-			message("File ", intensityFile, " does not exist.")

-			message("Not loading it, but running SNPRMA from scratch.")

-		}

 	if (!load.it){

 		res <- snprma(filenames,

 			      fitMixture=TRUE,

@@ -24,78 +28,96 @@ crlmm <- function(filenames, row.names=TRUE, col.names=TRUE,

 			      sns=sns,

 			      AFile=AFile,

 			      BFile=BFile)

-		save(res, file=intensityFile)##file.path(dirname(intensityFile), "res.rda"))

+		save(res, file=snprmaFile)##file.path(dirname(snprmaFile), "res.rda"))

 	} else {

-		message("Loading quantile normalized intensities...")

-		load(intensityFile)

+		message("Loading snprmaFile...")

+		load(snprmaFile)

 		res <- get("res")

 	}

+	SKW <- res[["SKW"]]

+	SNR <- res[["SNR"]]	

 	load(AFile)

 	load(BFile)

 	if(isPackageLoaded("ff")) {open(A); open(B)}

 	if(row.names) row.names <- res$gns

 	if(col.names) col.names <- res$sns	

 ##	}else{

-##		if (verbose) message("Loading ", intensityFile, ".")

-##		obj <- load(intensityFile)

+##		if (verbose) message("Loading ", snprmaFile, ".")

+##		obj <- load(snprmaFile)

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

 ##		if (obj != "res")

-##			stop("Object in ", intensityFile, " seems to be invalid.")

+##			stop("Object in ", snprmaFile, " seems to be invalid.")

 ##	}

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

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

 	gns <- res[["gns"]]

 	sns <- colnames(A)

+	nc <- ncol(A)	

 	##

 	##Ensures that if ff package is loaded, ordinary matrices are still passed to crlmmGT

-	aMatrix <- A[1:length(gns), ]

-	bMatrix <- B[1:length(gns), ]

-	if(isPackageLoaded("ff")) close(B)

-	rm(B); gc()

-	##	res2 <- crlmmGT(res[["A"]], res[["B"]], res[["SNR"]],

-	res2 <- crlmmGT(aMatrix, bMatrix, res[["SNR"]],

-			res[["mixtureParams"]], res[["cdfName"]],

-			gender=gender, row.names=row.names,

-			col.names=col.names, recallMin=recallMin,

-			recallRegMin=1000, SNRMin=SNRMin,

-			returnParams=returnParams, badSNP=badSNP,

-			verbose=verbose)

-	rm(aMatrix, bMatrix); gc()

-	res2[["SNR"]] <- res[["SNR"]]

-	res2[["SKW"]] <- res[["SKW"]]

+	if(nc > BS){

+		## genotype the samples in batches

+

+		##number crlmm batches

+		N <- ceiling(nc/BS)

+		##samples per batch

+		S <- ceiling(nc/N)

+		colindex <- split(1:nc, rep(1:nc, each=S, length.out=nc))

+	} else {

+		colindex <- list(1:nc)

+	}

+	if(length(colindex) > 1)

+		message("Calling genotypes in batches of size ", length(colindex[[1]]), " to reduce required RAM")	

 	if(isPackageLoaded("ff")){

-		if(file.exists(callsFile)) unlink(callsFile)

-		calls <- ff(dim=c(nrow(A), ncol(A)),

-			    vmode="integer",

-			    finalizer="close",

-			    dimnames=dimnames(A),

-			    overwrite=TRUE)

+		confs <- initializeBigMatrix(nrow(A), ncol(A))

+		calls <- initializeBigMatrix(nrow(A), ncol(A))

 	} else{

-		calls <- matrix(NA, nrow(A), ncol(A), dimnames=dimnames(A))

+		confs <- matrix(NA, nrow(A), ncol(A))

+		calls <- matrix(NA, nrow(A), ncol(A))

+	}	

+	sex <- vector("list", length(colindex))

+	for(i in seq(along=colindex)){

+		aMatrix <- A[1:length(gns), colindex[[i]]]

+		bMatrix <- B[1:length(gns), colindex[[i]]]

+		snr <- res[["SNR"]][colindex[[i]]]

+		mix <- res[["mixtureParameters"]][, colindex[[i]]]

+		##	res2 <- crlmmGT(res[["A"]], res[["B"]], res[["SNR"]],

+		res2 <- tryCatch(crlmmGT(aMatrix,

+					 bMatrix,

+					 snr,

+					 mix,

+					 res[["cdfName"]],

+					 gender=gender,

+					 row.names=row.names,

+					 col.names=col.names[colindex[[i]]],

+					 recallMin=recallMin,

+					 recallRegMin=1000,

+					 SNRMin=SNRMin,

+					 returnParams=returnParams,

+					 badSNP=badSNP,

+					 verbose=verbose), error=function(e) NULL)

+		if(is.null(res2)) {

+			unlink(callsFile)

+			unlink(confsFile)

+			stop("error in call to crlmmGT")

+		}

+		rm(aMatrix, bMatrix, snr, mix); gc()

+		calls[1:length(gns), colindex[[i]]] <- res2[["calls"]]

+		confs[1:length(gns), colindex[[i]]] <- res2[["confs"]]

+		sex[[i]] <- res2[["gender"]]

+		rm(res2); gc()

 	}

-	calls[1:length(gns), ] <- res2[["calls"]]

+	if(isPackageLoaded("ff")) close(B)

+	rm(B); gc()

 	if(isPackageLoaded("ff")) close(calls)

 	save(calls, file=callsFile)

 	rm(calls); gc()

-	if(isPackageLoaded("ff")){

-		if(file.exists(confsFile)) unlink(confsFile)		

-		confs <- ff(dim=c(nrow(A), ncol(A)),

-			    vmode="integer",

-			    finalizer="close",

-			    dimnames=dimnames(A),

-			    overwrite=TRUE)

-	} else{

-		confs <- matrix(NA, nrow(A), ncol(A), dimnames=dimnames(A))

-	}

-	confs[1:length(gns), ] <- res2[["confs"]]

 	if(isPackageLoaded("ff")) close(confs)

 	save(confs, file=confsFile)

 	rm(confs); gc()

-	sampleStats <- data.frame(SNR=res2[["SNR"]],

-				  SKW=res2[["SKW"]],

-				  gender=res2[["gender"]])

-	save(sampleStats, file=file.path(dirname(confsFile), "sampleStats.rda"))

-	if(isPackageLoaded("ff")){ close(A)}

+	res$gender <- unlist(sex)

+	save(res, file=snprmaFile)

+	if(isPackageLoaded("ff")) close(A)

 	rm(list=ls()); gc()

 	return()

 }

@@ -50,13 +50,13 @@ snprma <- function(filenames, mixtureSampleSize=10^5,

 	cnProbes <- get("cnProbes")

 	nr <- length(pnsa) + nrow(cnProbes)

 	nc <- length(filenames)

-	dns <- list(c(gns, rownames(cnProbes)), sns)  

+	dns <- list(c(gns, rownames(cnProbes)), sns)

 	if(isPackageLoaded("ff")){

 		message("ff package is loaded.  Creating ff objects for storing normalized and summarized intensities")

-		A <- ff(dim=c(nr, nc), vmode="integer", finalizer="close", dimnames=dns,

-			overwrite=TRUE)

-		B <- ff(dim=c(nr, nc), vmode="integer", finalizer="close", dimnames=dns,

-			overwrite=TRUE)

+		A <- initializeBigMatrix(nr, nc)

+		B <- initializeBigMatrix(nr, nc)

+		open(A)

+		open(B)

 	}  else {

 		A <- matrix(0L, nr, nc, dimnames=dns)

 		B <- matrix(0L, nr, nc, dimnames=dns)

@@ -98,7 +98,8 @@ snprma <- function(filenames, mixtureSampleSize=10^5,

 	if (!fitMixture) SNR <- mixtureParams <- NA

 	save(A, file=AFile)

 	save(B, file=BFile)

-	return(list(sns=sns, gns=gns, SNR=SNR, SKW=SKW, mixtureParams=mixtureParams, cdfName=cdfName))

+	return(list(sns=sns, gns=gns, SNR=SNR, SKW=SKW, mixtureParams=mixtureParams, cdfName=cdfName,

+		    cnnames=rownames(cnProbes)))

 }

 fitAffySnpMixture56 <- function(S, M, knots, probs=rep(1/3, 3), eps=.01, maxit=10, verbose=FALSE){

@@ -183,4 +183,39 @@ isValidCdfName <- function(cdfName){

 	return(result)

 }

+initializeBigMatrix <- function(nr, nc, batch, vmode="integer"){

+	if(prod(nr, nc) > 2^31){

+		##Need multiple matrices

+		## -- use ffdf

+			

+		## How many samples per ff object

+		S <- floor(2^31/nr - 1)

+

+		## How many ff objects

+		L <- ceiling(nc/S)

+

+		results <- vector("list", L)

+		##resultsB <- vector("list", L)			

+		for(i in 1:(L-1)){  ## the Lth object may have fewer than nc columns

+			results[[i]] <- ff(dim=c(nr, S),

+					   vmode=vmode,

+					   finalizer="close",

+					   overwrite=TRUE)

+		}

+		##the Lth element

+		leftOver <- nc - ((L-1)*S)

+		results[[L]] <- ff(dim=c(nr, leftOver),

+				   vmode=vmode,

+				   finalizer="close",

+				   overwrite=TRUE)

+		resultsff <- do.call(ffdf, results)

+	} else {

+		resultsff <- ff(dim=c(nr, nc),

+				vmode=vmode,

+				finalizer="close",

+				overwrite=TRUE)

+	}

+	return(resultsff)

+}

+

deleted file mode 100644

@@ -1,24 +0,0 @@

-\name{AlleleSet-methods}

-\docType{methods}

-

-\title{Indicator for polymorphic probes}

-\description{

-  

-  This functions uses the annotation slot of the AlleleSet object to load

-  the corresponding annotation package and determine whether each probe

-  in the object interrogates a polymorphic or nonpolymorphic allele. For

-  instance, in the Affy 6.0 platform roughly 900,000 of the 1.8 million

-  markers are for polymorphic alleles.

-  

-}

-\usage{

-A(object)

-B(object)

-}

-\arguments{

-  \item{object}{AlleleSet object}

-}

-\value{

-	matrix of normalized intensities

-}

-\keyword{manip}

@@ -8,14 +8,16 @@

 }

 \usage{

     cnOptions(outdir = "./", cdfName, 

-              AFile="A",

-	      BFile="B",

-	      callsFile="genotypes",

+              AFile="A.rda",

+	      BFile="B.rda",

+	      CAFile="CA.rda",

+	      CBFile="CB.rda",

+	      callsFile="genotypes.rda",

 	      rgFile="rgFile.rda",

 	      cnFile="cnSet",

 	      cnrmaFile="cn_rmaResult.rda",

 	      snprmaFile="snp_rmaResult.rda",

-	      confsFile="confsFile",

+	      confsFile="confsFile.rda",

 	      save.it = TRUE, 

 	      save.cnset = TRUE, 

 	      load.it = TRUE, 

@@ -25,7 +27,8 @@

 	      GT.CONF.THR = 0.99, PHI.THR = 2^6, nHOM.THR = 5, 

 	      MIN.NU = 2^3, MIN.PHI = 2^3, 

 	      THR.NU.PHI = TRUE, thresholdCopynumber = TRUE, unlink = TRUE, 

-xo	      use.poe=FALSE, crlmmBatchSize=1100, ...)

+	      use.poe=FALSE, 

+	      crlmmBatchSize=1000, ...)

 }

 \arguments{

   \item{outdir}{

@@ -45,6 +48,21 @@ xo	      use.poe=FALSE, crlmmBatchSize=1100, ...)

   File name for quantile normalized intensities for the B allele is <outdir>/BFile<.rda>

+}

+

+  \item{CAFile}{

+

+  File name for A allele copy number is <outdir>/AFile<.rda>.  Copy

+  number at nonpolymorphic loci is stored in this file.

+

+}

+

+  \item{CBFile}{

+

+  File name for B allele copy number is <outdir>/AFile<.rda>.  The 'B

+  allele' copy number at nonpolymorphic loci is zero (copy number for

+  nonpolymorphic loci is stored in the 'A allele').