Browse code

A number of changes to fit.lm1 and crlmmCopynumber. The number of genotypes is now stored as part of the CNSet class object.

git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/crlmm@48958 bc3139a8-67e5-0310-9ffc-ced21a209358

Rob Scharp authored on 21/08/2010 02:49:47
Showing8 changed files

... ...
@@ -29,7 +29,6 @@ Collate: AllGenerics.R
29 29
 	 AllClasses.R
30 30
 	 methods-CNSet.R
31 31
 	 methods-eSet.R
32
-	 methods-LinearModelParameter.R
33 32
          methods-SnpSuperSet.R
34 33
          cnrma-functions.R
35 34
          crlmm-functions.R
... ...
@@ -29,7 +29,7 @@ importClassesFrom(oligoClasses, SnpSuperSet, AlleleSet, CNSet)
29 29
 importMethodsFrom(oligoClasses, allele, calls, "calls<-", confs,
30 30
 		  "confs<-", cnConfidence, "cnConfidence<-", isSnp,
31 31
 		  chromosome, position, A, B,
32
-		  "A<-", "B<-", open, close, lM, "lM<-", flags, "flags<-")
32
+		  "A<-", "B<-", open, close, lM, "lM<-", flags)
33 33
 
34 34
 importFrom(oligoClasses, chromosome2integer, celfileDate, list.celfiles,
35 35
            copyNumber, initializeBigMatrix, initializeBigVector)
... ...
@@ -63,19 +63,23 @@ export(crlmm,
63 63
        crlmmIllumina,
64 64
        crlmmIllumina2,
65 65
        ellipseCenters,
66
+       genotype,
66 67
        readIdatFiles,
67 68
        readIdatFiles2,
68 69
        snprma,
69 70
        snprma2,
70 71
        crlmm2,
71 72
        genotype2, genotypeLD,
72
-       crlmmCopynumber2, crlmmCopynumberLD)
73
+       crlmmCopynumber2, crlmmCopynumberLD, crlmmCopynumber)
73 74
 export(constructIlluminaCNSet)
74 75
 ##export(linesCNSet)
75
-export(computeCN, fit.lm1, fit.lm2, fit.lm3, fit.lm4, construct,
76
+export(fit.lm1, fit.lm2, fit.lm3, fit.lm4, construct,
76 77
        dqrlsWrapper, fit.wls)
77 78
 export(computeCopynumber, ACN)
78 79
 export(totalCopyNumber)
80
+export(cnrma, cnrma2)
81
+
82
+exportMethods("A<-", "B<-", "nuA", "nuB", "phiA", "phiB", "nuA<-", "nuB<-", "phiA<-", "phiB<-", "numberGenotype<-")
79 83
 
80 84
 ## For debugging
81 85
 ## exportPattern("^[^\\.]")
... ...
@@ -36,5 +36,8 @@ setGeneric("tau2B<-", function(object, value) standardGeneric("tau2B<-"))
36 36
 setGeneric("corrAA<-", function(object, value) standardGeneric("corrAA<-"))
37 37
 setGeneric("corrAB<-", function(object, value) standardGeneric("corrAB<-"))
38 38
 setGeneric("corrBB<-", function(object, value) standardGeneric("corrBB<-"))
39
+setGeneric("flags<-", function(object, value) standardGeneric("flags<-"))
40
+
41
+setGeneric("numberGenotype<-", function(object, value, ...) standardGeneric("numberGenotype<-"))
39 42
 
40 43
 
... ...
@@ -80,8 +80,9 @@ construct <- function(filenames, cdfName, copynumber=FALSE,
80 80
 		protocolData <- annotatedDataFrameFrom(A(cnSet), byrow=FALSE)
81 81
 	}
82 82
 	rownames(pData(protocolData)) <- sns
83
-	protocolData(callSet) <- protocolData
83
+	protocolData(cnSet) <- protocolData
84 84
 	featureData(cnSet) <- featureData
85
+	featureNames(cnSet) <- featureNames(featureData)
85 86
 	pd <- data.frame(matrix(NA, nc, 3), row.names=sns)
86 87
 	colnames(pd)=c("SKW", "SNR", "gender")
87 88
 	phenoData(cnSet) <- new("AnnotatedDataFrame", data=pd)
... ...
@@ -179,7 +180,7 @@ construct <- function(filenames, cdfName, copynumber=FALSE,
179 180
 
180 181
 
181 182
 ##genotypeLargeData
182
-genotypeLD <- function(filenames,
183
+genotype <- function(filenames,
183 184
 		       cdfName,
184 185
 		       batch,
185 186
 		       mixtureSampleSize=10^5,
... ...
@@ -196,109 +197,125 @@ genotypeLD <- function(filenames,
196 197
 		       gender=NULL,
197 198
 		       returnParams=TRUE,
198 199
 		       badSNP=0.7){
199
-	if(!isPackageLoaded("ff")) stop("Must load package 'ff'")
200
+	is.lds <- ifelse(isPackageLoaded("ff"), TRUE, FALSE)
200 201
 	if(!copynumber){
201
-		callSet <- crlmm2(filenames=filenames,
202
-				  cdfName=cdfName,
203
-				  mixtureSampleSize=mixtureSampleSize,
204
-				  eps=eps,
205
-				  verbose=verbose,
206
-				  sns=sns,
207
-				  probs=probs,
208
-				  DF=DF,
209
-				  SNRMin=SNRMin,
210
-				  recallMin=recallMin,
211
-				  recallRegMin=recallRegMin,
212
-				  gender=gender,
213
-				  returnParams=returnParams,
214
-				  badSNP=badSNP)
202
+		FUN <- ifelse(is.lds, "crlmm2", "crlmm")
203
+		callSet <- FUN(filenames=filenames,
204
+			       cdfName=cdfName,
205
+			       mixtureSampleSize=mixtureSampleSize,
206
+			       eps=eps,
207
+			       verbose=verbose,
208
+			       sns=sns,
209
+			       probs=probs,
210
+			       DF=DF,
211
+			       SNRMin=SNRMin,
212
+			       recallMin=recallMin,
213
+			       recallRegMin=recallRegMin,
214
+			       gender=gender,
215
+			       returnParams=returnParams,
216
+			       badSNP=badSNP)
215 217
 		return(callSet)
216 218
 	}
217 219
 	if(missing(cdfName)) stop("must specify cdfName")
218 220
 	if(!isValidCdfName(cdfName)) stop("cdfName not valid.  see validCdfNames")
219 221
 	if(missing(sns)) sns <- basename(filenames)
220
-	## callSet contains potentially very big matrices
221 222
 	callSet <- construct(filenames=filenames,
222 223
 			     cdfName=cdfName,
223 224
 			     copynumber=TRUE,
224 225
 			     sns=sns,
225 226
 			     verbose=verbose,
226 227
 			     batch=batch)
227
-
228
+	open(callSet)
228 229
 	mixtureParams <- matrix(NA, 4, length(filenames))
229
-	snp.index <- which(isSnp(callSet)==1)
230
-	snprmaRes <- snprma2(filenames=filenames,
231
-			     mixtureSampleSize=mixtureSampleSize,
232
-			     fitMixture=TRUE,
233
-			     eps=eps,
234
-			     verbose=verbose,
235
-			     seed=seed,
236
-			     cdfName=cdfName,
237
-			     sns=sns)
230
+	is.snp <- isSnp(callSet)
231
+	snp.index <- which(is.snp)
232
+	FUN <- ifelse(is.lds, "snprma2", "snprma")
233
+	snprmaFxn <- function(FUN,...){
234
+		switch(FUN,
235
+		       snprma=snprma(...),
236
+		       snprma2=snprma2(...))
237
+	}
238
+	snprmaRes <- snprmaFxn(FUN, filenames=filenames,
239
+			       mixtureSampleSize=mixtureSampleSize,
240
+			       fitMixture=TRUE,
241
+			       eps=eps,
242
+			       verbose=verbose,
243
+			       seed=seed,
244
+			       cdfName=cdfName,
245
+			       sns=sns)
238 246
 	if(verbose) message("Finished preprocessing.")
239
-	open(snprmaRes[["A"]])
240
-	open(snprmaRes[["B"]])
241
-	open(snprmaRes[["SNR"]])
242
-	open(snprmaRes[["SKW"]])
243
-	open(snprmaRes[["mixtureParams"]])
244
-	if(verbose) message("Updating elements of callSet")
245
-##	bb = ocProbesets()*ncol(A)*8
246
-##	ffrowapply(A(callSet)[i1:i2, ] <- snprmaRes[["A"]][i1:i2, ], X=snprmaRes[["A"]], BATCHBYTES=bb)
247
-##	ffrowapply(B(callSet)[i1:i2, ] <- snprmaRes[["B"]][i1:i2, ], X=snprmaRes[["B"]], BATCHBYTES=bb)
248
-	##batches <- splitIndicesByLength(1:nrow(snprmaRes[["A"]]), ocProbesets())
249
-	AA <- A(callSet)
250
-	BB <- B(callSet)
251
-	for(j in 1:ncol(callSet)){
252
-		AA[snp.index, j] <- snprmaRes[["A"]][, j]
253
-		BB[snp.index, j] <- snprmaRes[["B"]][, j]
247
+	if(is.lds){
248
+		open(snprmaRes[["A"]])
249
+		open(snprmaRes[["B"]])
250
+		bb = ocProbesets()*ncol(A)*8
251
+		ffrowapply(A(callSet)[i1:i2, ] <- snprmaRes[["A"]][i1:i2, ], X=snprmaRes[["A"]], BATCHBYTES=bb)
252
+		ffrowapply(B(callSet)[i1:i2, ] <- snprmaRes[["B"]][i1:i2, ], X=snprmaRes[["B"]], BATCHBYTES=bb)
253
+	} else{
254
+		A(callSet)[snp.index, ] <- snprmaRes[["A"]]
255
+		B(callSet)[snp.index, ] <- snprmaRes[["B"]]
254 256
 	}
255
-	if(verbose) message("Finished updating elements of callSet")
256
-	stopifnot(identical(featureNames(callSet)[snp.index], snprmaRes$gns))
257
+##	stopifnot(identical(featureNames(callSet)[snp.index], snprmaRes$gns))
257 258
 	pData(callSet)$SKW <- snprmaRes$SKW
258 259
 	pData(callSet)$SNR <- snprmaRes$SNR
259 260
 	mixtureParams <- snprmaRes$mixtureParams
260
-	np.index <- which(isSnp(callSet) == 0)
261
-	cnrmaRes <- cnrma2(A=AA,
262
-			   filenames=filenames,
263
-			   row.names=featureNames(callSet)[np.index],
264
-			   cdfName=cdfName,
265
-			   sns=sns,
266
-			   seed=seed,
267
-			   verbose=verbose)
268
-	rm(cnrmaRes); ##gc()
269
-	## as.matrix needed when ffdf is used
270
-	tmp <- crlmmGT2(A=snprmaRes[["A"]],
271
-			B=snprmaRes[["B"]],
272
-			SNR=snprmaRes[["SNR"]],
273
-			mixtureParams=snprmaRes[["mixtureParams"]],
274
-			cdfName=cdfName,
275
-			row.names=NULL, ##featureNames(callSet),##[snp.index],
276
-			col.names=sampleNames(callSet),
277
-			probs=probs,
278
-			DF=DF,
279
-			SNRMin=SNRMin,
280
-			recallMin=recallMin,
281
-			recallRegMin=recallRegMin,
282
-			gender=gender,
283
-			verbose=verbose,
284
-			returnParams=returnParams,
285
-			badSNP=badSNP)
286
-	if(verbose) message("Leaving crlmmGT2")
287
-	open(tmp[["calls"]])
288
-	open(tmp[["confs"]])
289
-	##bb = ocProbesets()*ncol(A)*8
290
-	GT <- snpCall(callSet)
291
-	GTP <- snpCallProbability(callSet)
292
-	for(j in 1:ncol(callSet)){
293
-		GT[snp.index, j] <- tmp[["calls"]][, j]
294
-		GTP[snp.index, j] <- tmp[["confs"]][, j]
261
+	np.index <- which(!is.snp)
262
+	if(verbose) message("Normalizing nonpolymorphic markers")
263
+	FUN <- ifelse(is.lds, "cnrma2", "cnrma")
264
+	## main purpose is to update 'alleleA'
265
+	cnrmaFxn <- function(FUN,...){
266
+		switch(FUN,
267
+		       cnrma=cnrma(...),
268
+		       cnrma2=cnrma2(...))
269
+	}
270
+	## consider passing only A for NPs.
271
+	AA <- cnrmaFxn(FUN, A=A(callSet),
272
+		       filenames=filenames,
273
+		       row.names=featureNames(callSet)[np.index],
274
+		       cdfName=cdfName,
275
+		       sns=sns,
276
+		       seed=seed,
277
+		       verbose=verbose)
278
+	if(!is.lds) A(callSet) <- AA
279
+	rm(AA)
280
+	FUN <- ifelse(is.lds, "crlmmGT2", "crlmmGT")
281
+	## genotyping
282
+	crlmmGTfxn <- function(FUN,...){
283
+		switch(FUN,
284
+		       crlmmGT2=crlmmGT2(...),
285
+		       crlmmGT=crlmmGT(...))
286
+	}
287
+	tmp <- crlmmGTfxn(FUN,
288
+			  A=snprmaRes[["A"]],
289
+			  B=snprmaRes[["B"]],
290
+			  SNR=snprmaRes[["SNR"]],
291
+			  mixtureParams=snprmaRes[["mixtureParams"]],
292
+			  cdfName=cdfName,
293
+			  row.names=NULL,
294
+			  col.names=sampleNames(callSet),
295
+			  probs=probs,
296
+			  DF=DF,
297
+			  SNRMin=SNRMin,
298
+			  recallMin=recallMin,
299
+			  recallRegMin=recallRegMin,
300
+			  gender=gender,
301
+			  verbose=verbose,
302
+			  returnParams=returnParams,
303
+			  badSNP=badSNP)
304
+	if(verbose) message("Genotyping finished.  Updating container with genotype calls and confidence scores.")
305
+	if(is.lds){
306
+		open(tmp[["calls"]])
307
+		open(tmp[["confs"]])
308
+		ffrowapply(snpCall(callSet)[i1:i2, ] <- tmp[["calls"]][i1:i2, ], X=tmp[["calls"]], BATCHBYTES=bb)
309
+		ffrowapply(snpCallProbability(callSet)[i1:i2, ] <- tmp[["confs"]][i1:i2, ], X=tmp[["confs"]], BATCHBYTES=bb)
310
+	} else {
311
+		calls(callSet)[snp.index, ] <- tmp[["calls"]]
312
+		snpCallProbability(callSet)[snp.index, ] <- tmp[["confs"]]
295 313
 	}
296
-##	ffrowapply(snpCall(callSet)[i1:i2, ] <- tmp[["calls"]][i1:i2, ], X=tmp[["calls"]], BATCHBYTES=bb)
297
-##	ffrowapply(snpCallProbability(callSet)[i1:i2, ] <- tmp[["confs"]][i1:i2, ], X=tmp[["confs"]], BATCHBYTES=bb)
298 314
 	callSet$gender <- tmp$gender
315
+	close(callSet)
299 316
 	return(callSet)
300 317
 }
301
-genotype <- genotype2 <- genotypeLD
318
+genotypeLD <- genotype2 <- genotype
302 319
 
303 320
 rowCovs <- function(x, y, ...){
304 321
 	notna <- !is.na(x)
... ...
@@ -763,11 +780,10 @@ crlmmCopynumberLD <- function(object,
763 780
 }
764 781
 crlmmCopynumber2 <- crlmmCopynumberLD
765 782
 
766
-fit.lm1 <- function(idxBatch,
767
-		    snpBatches,
783
+fit.lm1 <- function(strata.index,
784
+		    index.list,
768 785
 		    marker.index,
769 786
 		    object,
770
-		    Ns,
771 787
 		    batchSize,
772 788
 		    SNRMin,
773 789
 		    MIN.SAMPLES,
... ...
@@ -779,25 +795,26 @@ fit.lm1 <- function(idxBatch,
779 795
 		    MIN.PHI,
780 796
 		    verbose,...){
781 797
 	if(isPackageLoaded("ff")) physical <- get("physical")
782
-	if(verbose) message("Probe stratum ", idxBatch, " of ", length(snpBatches))
783
-	snps <- snpBatches[[idxBatch]]
784
-	batches <- split(seq(along=batch(object)), batch(object))
798
+	is.lds <- ifelse(is(calls(object), "ffdf") | is(calls(object), "ff_matrix"), TRUE, FALSE)
799
+	if(verbose) message("Probe stratum ", strata.index, " of ", length(index.list))
800
+	snps <- index.list[[strata.index]]
801
+	batches <- split(seq_along(batch(object)), batch(object))
785 802
 	batches <- batches[sapply(batches, length) >= MIN.SAMPLES]
786
-
787 803
 	open(object)
788
-	open(snpflags)
789
-##	open(normal)
790
-
791 804
 	corrAB <- corrBB <- corrAA <- sig2B <- sig2A <- tau2B <- tau2A <- matrix(NA, length(snps), length(unique(batch(object))))
792 805
 	flags <- nuA <- nuB <- phiA <- phiB <- corrAB
806
+	## should the 'normal' indicator
807
+	NORM <- 1
793 808
 
794
-##	normal.snps <- normal[snps, ]
809
+	Ns.names <- ls(numberGenotype(object))
810
+	batchnames <- batchNames(object)
795 811
 	GG <- as.matrix(calls(object)[snps, ])
796 812
 	CP <- as.matrix(snpCallProbability(object)[snps, ])
797 813
 	AA <- as.matrix(A(object)[snps, ])
798 814
 	BB <- as.matrix(B(object)[snps, ])
799 815
 	zzzz <- 1
800 816
 	for(k in batches){
817
+		this.batch <- unique(as.character(batch(object)[k]))
801 818
 		zzzz <- zzzz+1
802 819
 		G <- GG[, k]
803 820
 		##NORM <- normal.snps[, k]
... ...
@@ -816,7 +833,7 @@ fit.lm1 <- function(idxBatch,
816 833
 		vA <- shrink(vA, Ns, DF.PRIOR)
817 834
 		vB <- shrink(vB, Ns, DF.PRIOR)
818 835
 		##location and scale
819
-		J <- match(unique(batch(object)[k]), unique(batch(object)))
836
+		J <- match(unique(batch(object)[k]), batchnames)##unique(batch(object)))
820 837
 		##background variance for alleleA
821 838
 		taus <- applyByGenotype(log2(A), rowMAD, G)^2
822 839
 		tau2A[, J] <- shrink(taus[, 3, drop=FALSE], Ns[, 3], DF.PRIOR)
... ...
@@ -877,7 +894,6 @@ fit.lm1 <- function(idxBatch,
877 894
 		negB <- rowSums(muB < 0) > 0
878 895
 		flags[, J] <- rowSums(Ns == 0) > 0
879 896
 		##flags[, J] <- index[[1]] | index[[2]] | index[[3]] | rowSums(
880
-
881 897
 		## we're regressing on the medians using the standard errors (hence the division by N) as weights
882 898
 		##formerly coefs()
883 899
 		Np <- Ns
... ...
@@ -904,17 +920,21 @@ fit.lm1 <- function(idxBatch,
904 920
 		nuB[, J] <- res[1, ]
905 921
 		##phiB[, J] <- res[[2]]
906 922
 		phiB[, J] <- res[2, ]
907
-		if(THR.NU.PHI){
908
-			nuA[nuA[, J] < MIN.NU, J] <- MIN.NU
909
-			nuB[nuB[, J] < MIN.NU, J] <- MIN.NU
910
-			phiA[phiA[, J] < MIN.PHI, J] <- MIN.PHI
911
-			phiB[phiB[, J] < MIN.PHI, J] <- MIN.PHI
912
-		}
913 923
 ##		cA[, k] <- matrix((1/phiA[, J]*(A-nuA[, J])), nrow(A), ncol(A))
914 924
 ##		cB[, k] <- matrix((1/phiB[, J]*(B-nuB[, J])), nrow(B), ncol(B))
915
-		rm(G, A, B, NORM, wA, wB, YA,YB, res, negA, negB, Np, Ns)
925
+		jj <- match(this.batch, Ns.names)
926
+		numberGenotype(object, this.batch)[snps, ] <- Ns
927
+		rm(G, A, B, wA, wB, YA,YB, res, negA, negB, Np, Ns)
916 928
 		##gc()
917 929
 	}
930
+	nGt <- lapply(nGt, function(x){ colnames(x) <- c("AA", "AB", "BB"); return(x)})
931
+
932
+	if(THR.NU.PHI){
933
+		nuA[nuA < MIN.NU] <- MIN.NU
934
+		nuB[nuB < MIN.NU] <- MIN.NU
935
+		phiA[phiA < MIN.PHI] <- MIN.PHI
936
+		phiB[phiB < MIN.PHI] <- MIN.PHI
937
+	}
918 938
 ##	cA[cA < 0.05] <- 0.05
919 939
 ##	cB[cB < 0.05] <- 0.05
920 940
 ##	cA[cA > 5] <-  5
... ...
@@ -923,54 +943,30 @@ fit.lm1 <- function(idxBatch,
923 943
 ##	cB <- matrix(as.integer(cB*100), nrow(cB), ncol(cB))
924 944
 ##	CA(object)[snps, ] <- cA
925 945
 ##	CB(object)[snps, ] <- cB
926
-
927
-	snpflags[snps, ] <- flags
928
-	lapply(lM(object), open)
929
-
930
-	tmp <- physical(lM(object))$tau2A
931
-	tmp[snps, ] <- tau2A
932
-	lM(object)$tau2A <- tmp
933
-	tmp <- physical(lM(object))$tau2B
934
-	tmp[snps, ] <- tau2B
935
-	lM(object)$tau2B <- tmp
936
-	tmp <- physical(lM(object))$tau2B
937
-	tmp[snps, ] <- tau2B
938
-	lM(object)$tau2B <- tmp
939
-	tmp <- physical(lM(object))$sig2A
940
-	tmp[snps, ] <- sig2A
941
-	lM(object)$sig2A <- tmp
942
-	tmp <- physical(lM(object))$sig2B
943
-	tmp[snps, ] <- sig2B
944
-	lM(object)$sig2B <- tmp
945
-	tmp <- physical(lM(object))$nuA
946
-	tmp[snps, ] <- nuA
947
-	lM(object)$nuA <- tmp
948
-	tmp <- physical(lM(object))$nuB
949
-	tmp[snps, ] <- nuB
950
-	lM(object)$nuB <- tmp
951
-	tmp <- physical(lM(object))$phiA
952
-	tmp[snps, ] <- phiA
953
-	lM(object)$phiA <- tmp
954
-	tmp <- physical(lM(object))$phiB
955
-	tmp[snps, ] <- phiB
956
-	lM(object)$phiB <- tmp
957
-	tmp <- physical(lM(object))$corrAB
958
-	tmp[snps, ] <- corrAB
959
-	lM(object)$corrAB <- tmp
960
-	tmp <- physical(lM(object))$corrAA
961
-	tmp[snps, ] <- corrAA
962
-	lM(object)$corrAA <- tmp
963
-	tmp <- physical(lM(object))$corrBB
964
-	tmp[snps, ] <- corrBB
965
-	lM(object)$corrBB <- tmp
966
-
967
-	lapply(assayData(object), close)
968
-	lapply(lM(object), close)
969
-	TRUE
946
+	if(is.lds) lapply(lM(object), open)
947
+	flags(object)[snps, ] <- flags
948
+	tau2A(object) <- tau2A
949
+	tau2B(object) <- tau2B
950
+	sigma2A(object) <- sig2A
951
+	sigma2B(object) <- sig2B
952
+	nuA(object) <- nuA
953
+	nuB(object) <- nuB
954
+	phiA(object) <- phiA
955
+	phiB(object) <- phiB
956
+	corrAA(object) <- corrAA
957
+	corrBB(object) <- corrBB
958
+	corrAB(object) <- corrAB
959
+	if(is.lds) {
960
+		lapply(assayData(object), close)
961
+		lapply(lM(object), close)
962
+		return(TRUE)
963
+	} else{
964
+		return(object)
965
+	}
970 966
 }
971 967
 
972
-fit.lm2 <- function(idxBatch,
973
-		    snpBatches,
968
+fit.lm2 <- function(strata.index,
969
+		    index.list,
974 970
 		    marker.index,
975 971
 		    object,
976 972
 		    Ns,
... ...
@@ -985,8 +981,8 @@ fit.lm2 <- function(idxBatch,
985 981
 		    MIN.PHI,
986 982
 		    verbose,...){
987 983
 	physical <- get("physical")
988
-	if(verbose) message("Probe stratum ", idxBatch, " of ", length(snpBatches))
989
-	snps <- snpBatches[[idxBatch]]
984
+	if(verbose) message("Probe stratum ", strata.index, " of ", length(index.list))
985
+	snps <- index.list[[strata.index]]
990 986
 	batches <- split(seq(along=batch(object)), batch(object))
991 987
 	batches <- batches[sapply(batches, length) >= MIN.SAMPLES]
992 988
 
... ...
@@ -1078,8 +1074,8 @@ fit.lm2 <- function(idxBatch,
1078 1074
 }
1079 1075
 
1080 1076
 
1081
-fit.lm3 <- function(idxBatch,
1082
-		    snpBatches,
1077
+fit.lm3 <- function(strata.index,
1078
+		    index.list,
1083 1079
 		    marker.index,
1084 1080
 		    object,
1085 1081
 		    Ns,
... ...
@@ -1094,8 +1090,8 @@ fit.lm3 <- function(idxBatch,
1094 1090
 		    MIN.PHI,
1095 1091
 		    verbose, ...){
1096 1092
 	physical <- get("physical")
1097
-	if(verbose) message("Probe stratum ", idxBatch, " of ", length(snpBatches))
1098
-		snps <- snpBatches[[idxBatch]]
1093
+	if(verbose) message("Probe stratum ", strata.index, " of ", length(index.list))
1094
+		snps <- index.list[[strata.index]]
1099 1095
 	batches <- split(seq(along=batch(object)), batch(object))
1100 1096
 	batches <- batches[sapply(batches, length) >= MIN.SAMPLES]
1101 1097
 
... ...
@@ -1315,8 +1311,8 @@ fit.lm3 <- function(idxBatch,
1315 1311
 	TRUE
1316 1312
 }
1317 1313
 
1318
-fit.lm4 <- function(idxBatch,
1319
-		    snpBatches,
1314
+fit.lm4 <- function(strata.index,
1315
+		    index.list,
1320 1316
 		    marker.index,
1321 1317
 		    object,
1322 1318
 		    Ns,
... ...
@@ -1331,11 +1327,11 @@ fit.lm4 <- function(idxBatch,
1331 1327
 		    MIN.PHI,
1332 1328
 		    verbose, ...){
1333 1329
 	physical <- get("physical")
1334
-	if(verbose) message("Probe stratum ", idxBatch, " of ", length(snpBatches))
1330
+	if(verbose) message("Probe stratum ", strata.index, " of ", length(index.list))
1335 1331
 	open(object)
1336 1332
 	open(normal)
1337 1333
 	open(snpflags)
1338
-	snps <- snpBatches[[idxBatch]]
1334
+	snps <- index.list[[strata.index]]
1339 1335
 	batches <- split(seq(along=batch(object)), batch(object))
1340 1336
 	batches <- batches[sapply(batches, length) >= MIN.SAMPLES]
1341 1337
 	nuA <- phiA <- sig2A <- tau2A <- matrix(NA, length(snps), length(unique(batch(object))))
... ...
@@ -1478,7 +1474,7 @@ whichPlatform <- function(cdfName){
1478 1474
 	return(platform)
1479 1475
 }
1480 1476
 
1481
-cnrma2 <- function(A, filenames, row.names, verbose=TRUE, seed=1, cdfName, sns){
1477
+cnrma <- function(A, filenames, row.names, verbose=TRUE, seed=1, cdfName, sns){
1482 1478
 	if(missing(cdfName)) stop("must specify cdfName")
1483 1479
 	pkgname <- getCrlmmAnnotationName(cdfName)
1484 1480
 	require(pkgname, character.only=TRUE) || stop("Package ", pkgname, " not available")
... ...
@@ -1486,24 +1482,48 @@ cnrma2 <- function(A, filenames, row.names, verbose=TRUE, seed=1, cdfName, sns){
1486 1482
 	if(verbose) message("Loading annotations for nonpolymorphic probes")
1487 1483
         loader("npProbesFid.rda", .crlmmPkgEnv, pkgname)
1488 1484
 	fid <- getVarInEnv("npProbesFid")
1489
-	SKW <- initializeBigVector("crlmmSKW.NP-", length(filenames), "double")
1490
-  	##A <- initializeBigMatrix("crlmmNP-", length(fid), length(filenames), "integer")
1485
+
1486
+	if(cdfName=="genomewidesnp6"){
1487
+		loader("1m_reference_cn.rda", .crlmmPkgEnv, pkgname)
1488
+	}
1489
+	if(cdfName=="genomewidesnp5"){
1490
+		loader("5.0_reference_cn.rda", .crlmmPkgEnv, pkgname)
1491
+	}
1492
+	reference <- getVarInEnv("reference")
1493
+	fid <- fid[match(row.names, names(fid))]
1494
+	np.index <- match(row.names, rownames(A))
1495
+	gns <- names(fid)
1496
+	set.seed(seed)
1497
+	idx2 <- sample(length(fid), 10^5)
1498
+	for (k in seq_along(filenames)){
1499
+		y <- as.matrix(read.celfile(filenames[k], intensity.means.only=TRUE)[["INTENSITY"]][["MEAN"]][fid])
1500
+		A[np.index, k] <- as.integer(normalize.quantiles.use.target(y, target=reference))
1501
+	}
1502
+	return(A)
1503
+}
1504
+
1505
+cnrma2 <- function(A, filenames, row.names, verbose=TRUE, seed=1, cdfName, sns){
1506
+	if(missing(cdfName)) stop("must specify cdfName")
1507
+	pkgname <- getCrlmmAnnotationName(cdfName)
1508
+	require(pkgname, character.only=TRUE) || stop("Package ", pkgname, " not available")
1509
+	if (missing(sns)) sns <- basename(filenames)
1491 1510
 	sampleBatches <- splitIndicesByNode(seq(along=filenames))
1492 1511
 	if(verbose) message("Processing nonpolymorphic probes for ", length(filenames), " files.")
1512
+	## updates A
1493 1513
 	ocLapply(sampleBatches,
1494 1514
 		 processCEL2,
1495 1515
 		 row.names=row.names,
1496 1516
 		 filenames=filenames,
1497 1517
 		 A=A,
1498
-		 SKW=SKW,
1499 1518
 		 seed=seed,
1500 1519
 		 pkgname=pkgname,
1501 1520
 		 cdfName=cdfName,
1502 1521
 		 neededPkgs=c("crlmm", pkgname))
1503
-	list(sns=sns, gns=row.names, SKW=SKW, cdfName=cdfName)
1522
+	##list(sns=sns, gns=row.names, SKW=SKW, cdfName=cdfName)
1523
+	return(A)
1504 1524
 }
1505 1525
 
1506
-processCEL2 <- function(i, filenames, row.names, A, SKW, seed, cdfName, pkgname){
1526
+processCEL2 <- function(i, filenames, row.names, A, seed, cdfName, pkgname){
1507 1527
 	if(cdfName=="genomewidesnp6"){
1508 1528
 		loader("1m_reference_cn.rda", .crlmmPkgEnv, pkgname)
1509 1529
 	}
... ...
@@ -1518,19 +1538,11 @@ processCEL2 <- function(i, filenames, row.names, A, SKW, seed, cdfName, pkgname)
1518 1538
 	gns <- names(fid)
1519 1539
 	set.seed(seed)
1520 1540
 	idx2 <- sample(length(fid), 10^5)
1521
-	open(A)
1522
-	open(SKW)
1523 1541
 	for (k in i){
1524 1542
 		y <- as.matrix(read.celfile(filenames[k], intensity.means.only=TRUE)[["INTENSITY"]][["MEAN"]][fid])
1525
-		x <- log2(y[idx2])
1526
-		SKW[k] <- mean((x-mean(x))^3)/(sd(x)^3)
1527
-		rm(x)
1528 1543
 		A[np.index, k] <- as.integer(normalize.quantiles.use.target(y, target=reference))
1529
-		rm(y)
1530 1544
 	}
1531
-	close(A)
1532
-	close(SKW)
1533
-	TRUE
1545
+	return(TRUE)
1534 1546
 }
1535 1547
 
1536 1548
 
... ...
@@ -2364,7 +2376,7 @@ imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
2364 2376
 ##}
2365 2377
 ##
2366 2378
 
2367
-##bias1 <- function(idxBatch,
2379
+##bias1 <- function(strata.index,
2368 2380
 ##		  snpBatches,
2369 2381
 ##		  index,
2370 2382
 ##		  object,
... ...
@@ -2376,7 +2388,7 @@ imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
2376 2388
 ##
2377 2389
 ##}
2378 2390
 
2379
-##bias2 <- function(idxBatch,
2391
+##bias2 <- function(strata.index,
2380 2392
 ##		  snpBatches,
2381 2393
 ##		  index,
2382 2394
 ##		  object,
... ...
@@ -2387,7 +2399,7 @@ imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
2387 2399
 ##	open(object)
2388 2400
 ##	open(normal)
2389 2401
 ##
2390
-##	nps <- snpBatches[[idxBatch]]
2402
+##	nps <- snpBatches[[strata.index]]
2391 2403
 ##	nuA <- lM(object)$nuA[nps, , drop=FALSE]
2392 2404
 ##	phiA <- lM(object)$phiA[nps, , drop=FALSE]
2393 2405
 ##	sig2A <- lM(object)$sig2A[nps, , drop=FALSE]
... ...
@@ -2767,91 +2779,95 @@ ellipseCenters <- function(object, index, allele, batch, log.it=TRUE){
2767 2779
 
2768 2780
 crlmmCopynumber <- function(object,
2769 2781
 			    filenames,
2770
-			    which.batches,
2771
-		      MIN.SAMPLES=10,
2772
-		      SNRMin=5,
2773
-		      MIN.OBS=1,
2774
-		      DF.PRIOR=50,
2775
-		      bias.adj=FALSE,
2776
-		      prior.prob=rep(1/4,4),
2777
-		      seed=1,
2778
-		      verbose=TRUE,
2779
-		      GT.CONF.THR=0.99,
2780
-		      PHI.THR=2^6,
2781
-		      nHOM.THR=5,
2782
-		      MIN.NU=2^3,
2783
-		      MIN.PHI=2^3,
2784
-		      THR.NU.PHI=TRUE,
2785
-		      thresholdCopynumber=TRUE,
2786
-		      type=c("autosome.snps", "autosome.nps", "X.snps", "X.nps")){
2787
-	stopifnot("batch" %in% varLabels(protocolData(object)))
2788
-	stopifnot("chromosome" %in% fvarLabels(object))
2789
-	stopifnot("position" %in% fvarLabels(object))
2790
-	stopifnot("isSnp" %in% fvarLabels(object))
2782
+			    MIN.SAMPLES=10,
2783
+			    SNRMin=5,
2784
+			    MIN.OBS=1,
2785
+			    DF.PRIOR=50,
2786
+			    bias.adj=FALSE,
2787
+			    prior.prob=rep(1/4,4),
2788
+			    seed=1,
2789
+			    verbose=TRUE,
2790
+			    GT.CONF.THR=0.99,
2791
+			    PHI.THR=2^6,
2792
+			    nHOM.THR=5,
2793
+			    MIN.NU=2^3,
2794
+			    MIN.PHI=2^3,
2795
+			    THR.NU.PHI=TRUE,
2796
+			    thresholdCopynumber=TRUE,
2797
+			    type=c("autosome.snps", "autosome.nps", "X.snps", "X.nps")){
2791 2798
 	batch <- batch(object)
2792 2799
 	is.snp <- isSnp(object)
2793 2800
 	is.autosome <- chromosome(object) < 23
2794 2801
 	is.annotated <- !is.na(chromosome(object))
2795 2802
 	is.X <- chromosome(object) == 23
2796
-	if(type=="autosome.snps"){
2797
-		marker.index <- which(is.snp & is.autosome & is.annotated)
2798
-		nr <- length(marker.index)
2799
-		Ns <- initializeBigMatrix("Ns", nr, 5)
2800
-		colnames(Ns) <- c("A", "B", "AA", "AB", "BB")
2801
-		if(verbose) message("Estimating allele-specific copy number at autosomal SNPs")
2802
-		FUN <- fit.lm1
2803
-	}
2804
-	if(type=="autosome.nps"){
2805
-		marker.index <- which(is.autosome & !is.snp & is.annotated)
2806
-		nr <- length(marker.index)
2807
-		Ns <- initializeBigMatrix("Ns", nr, 5)
2808
-		colnames(Ns) <- c("A", "B", "AA", "AB", "BB")
2809
-		if(verbose) message("Estimating allele-specific copy number at autosomal SNPs")
2810
-		FUN <- fit.lm2
2803
+	is.lds <- ifelse((is(calls(object), "ffdf") | is(calls(object), "ffdf")), TRUE, FALSE)
2804
+
2805
+	samplesPerBatch <- table(batch(object))
2806
+	if(any(samplesPerBatch < MIN.SAMPLES)){
2807
+		warning("The following batches have fewer than ", MIN.SAMPLES, ":")
2808
+		message(paste(samplesPerBatch[samplesPerBatch < MIN.SAMPLES], collapse=", "))
2809
+		message("Not estimating copy number for the above batches")
2811 2810
 	}
2812
-	if(type=="X.snps"){
2813
-		marker.index <- which(is.X & is.snp)
2814
-		nr <- length(marker.index)
2815
-		Ns <- initializeBigMatrix("Ns", nr, 5)
2816
-		colnames(Ns) <- c("A", "B", "AA", "AB", "BB")
2817
-		if(verbose) message("Estimating allele-specific copy number at autosomal SNPs")
2818
-		FUN <- fit.lm3
2811
+##	Ns <- initializeBigMatrix("Ns", nrow(object), 3)
2812
+##	colnames(Ns) <- c("AA", "AB", "BB")
2813
+	whichMarkers <- function(type, is.snp, is.autosome, is.annotated, is.X){
2814
+		switch(type,
2815
+		       autosome.snps=which(is.snp & is.autosome & is.annotated),
2816
+		       autosome.nps=which(!is.snp & is.autosome),
2817
+		       X.snps=which(is.snp & is.X),
2818
+		       X.nps=which(!is.snp & is.X),
2819
+		       stop("'type' must be one of 'autosome.snps', 'autosome.nps', 'X.snps', or 'X.nps'")
2820
+	       )
2819 2821
 	}
2820
-	if(type=="X.nps"){
2821
-		marker.index <- which(is.X & !is.snp)
2822
-		nr <- length(marker.index)
2823
-		Ns <- initializeBigMatrix("Ns", nr, 5)
2824
-		colnames(Ns) <- c("A", "B", "AA", "AB", "BB")
2825
-		if(verbose) message("Estimating allele-specific copy number at autosomal SNPs")
2826
-		FUN <- fit.lm4
2822
+	marker.index <- whichMarkers(type[[1]], is.snp, is.autosome, is.annotated, is.X)
2823
+	lmFxn <- function(type){
2824
+		switch(type,
2825
+		       autosome.snps="fit.lm1",
2826
+		       autosome.nps="fit.lm2",
2827
+		       X.snps="fit.lm3",
2828
+		       X.nps="fit.lm4")
2827 2829
 	}
2828
-	index.strata <- splitIndicesByLength(marker.index, ocProbesets())
2830
+	FUN <- lmFxn(type[[1]])
2831
+	index.list <- splitIndicesByLength(marker.index, ocProbesets())
2829 2832
 	## create a separate object for each strata and combine() at the very end.
2830
-	obj <- construct(filenames=filenames,
2831
-			 cdfName=annotation(object),
2832
-			 copynumber=TRUE,
2833
-			 batch=batch(object),
2834
-			 sns=sampleNames(object),
2835
-			 fns=featureNames(object)[marker.index])
2836
-	ocLapply(seq(along=index.strata),
2837
-		 FUN,
2838
-		 marker.index=marker.index,
2839
-		 object=obj,
2840
-		 Ns=Ns,
2841
-		 ##normal=normal,
2842
-		 ##snpflags=flags,
2843
-		 snpBatches=index.strata,
2844
-		 batchSize=ocProbesets(),
2845
-		 SNRMin=SNRMin,
2846
-		 MIN.SAMPLES=MIN.SAMPLES,
2847
-		 MIN.OBS=MIN.OBS,
2848
-		 DF=DF.PRIOR,
2849
-		 GT.CONF.THR=GT.CONF.THR,
2850
-		 THR.NU.PHI=THR.NU.PHI,
2851
-		 MIN.NU=MIN.NU,
2852
-		 MIN.PHI=MIN.PHI,
2853
-		 verbose=verbose,
2854
-		 neededPkgs="crlmm")
2833
+	## would need to put calls, confs, A, and B in each object.
2834
+	## can't assume that the dataset would be small enough to
2835
+	## represent as matrices.  would have to read from the
2836
+	## original matrix and reassign to a new matrix.  Probably not
2837
+	## worth the effort.
2838
+	if(is.lds){
2839
+		ocLapply(seq(along=index.list),
2840
+			 FUN,
2841
+			 index.list=index.list,
2842
+			 marker.index=marker.index,
2843
+			 object=object,
2844
+			 batchSize=ocProbesets(),
2845
+			 SNRMin=SNRMin,
2846
+			 MIN.SAMPLES=MIN.SAMPLES,
2847
+			 MIN.OBS=MIN.OBS,
2848
+			 DF=DF.PRIOR,
2849
+			 GT.CONF.THR=GT.CONF.THR,
2850
+			 THR.NU.PHI=THR.NU.PHI,
2851
+			 MIN.NU=MIN.NU,
2852
+			 MIN.PHI=MIN.PHI,
2853
+			 verbose=verbose,
2854
+			 neededPkgs="crlmm")
2855
+	} else {
2856
+		object <- fit.lm1(strata.index=1,
2857
+				  index.list=list(marker.index),
2858
+				  marker.index=marker.index,
2859
+				  object=object,
2860
+				  batchSize=batchSize,
2861
+				  SNRMin=SNRMin,
2862
+				  MIN.SAMPLES=MIN.SAMPLES,
2863
+				  MIN.OBS=MIN.OBS,
2864
+				  DF.PRIOR=DF.PRIOR,
2865
+				  GT.CONF.THR=GT.CONF.THR,
2866
+				  THR.NU.PHI=THR.NU.PHI,
2867
+				  MIN.NU=MIN.NU,
2868
+				  MIN.PHI=MIN.PHI,
2869
+				  verbose=verbose, ...)
2870
+	}
2855 2871
 	message("finished")
2856
-	return(obj)
2872
+	return(object)
2857 2873
 }
... ...
@@ -16,7 +16,6 @@ linearParamElementReplace <- function(obj, elt, value) {
16 16
     obj
17 17
 }
18 18
 
19
-
20 19
 setMethod("nuA", signature=signature(object="CNSet"), function(object) nu(object, "A"))
21 20
 setMethod("nuB", signature=signature(object="CNSet"), function(object) nu(object, "B"))
22 21
 setMethod("phiA", signature=signature(object="CNSet"), function(object) phi(object, "A"))
... ...
@@ -29,61 +28,76 @@ setMethod("corrAA", signature=signature(object="CNSet"), function(object) corr(o
29 28
 setMethod("corrAB", signature=signature(object="CNSet"), function(object) corr(object, "AB"))
30 29
 setMethod("corrBB", signature=signature(object="CNSet"), function(object) corr(object, "BB"))
31 30
 
32
-setReplaceMethod("nuA", signature=signature(object="CNSet", value="matrix"),
31
+setReplaceMethod("nuA", signature=signature(object="CNSet", value="ff_or_matrix"),
33 32
 	  function(object, value){
34 33
 		  linearParamElementReplace(object, "nuA", value)
35 34
 	  })
36 35
 
37
-setReplaceMethod("nuB", signature=signature(object="CNSet", value="matrix"),
36
+setReplaceMethod("nuB", signature=signature(object="CNSet", value="ff_or_matrix"),
38 37
 	  function(object, value){
39 38
 		  linearParamElementReplace(object, "nuB", value)
40 39
 })
41 40
 
42
-setReplaceMethod("phiA", signature=signature(object="CNSet", value="matrix"),
41
+setReplaceMethod("phiA", signature=signature(object="CNSet", value="ff_or_matrix"),
43 42
 	  function(object, value){
44 43
 		  linearParamElementReplace(object, "phiA", value)
45 44
 })
46 45
 
47
-setReplaceMethod("phiB", signature=signature(object="CNSet", value="matrix"),
46
+setReplaceMethod("phiB", signature=signature(object="CNSet", value="ff_or_matrix"),
48 47
 	  function(object, value){
49 48
 		  linearParamElementReplace(object, "phiB", value)
50 49
 })
51 50
 
52
-setReplaceMethod("sigma2A", signature=signature(object="CNSet", value="matrix"),
51
+setReplaceMethod("sigma2A", signature=signature(object="CNSet", value="ff_or_matrix"),
53 52
 	  function(object, value){
54 53
 		  linearParamElementReplace(object, "sig2A", value)
55 54
 })
56 55
 
57
-setReplaceMethod("sigma2B", signature=signature(object="CNSet", value="matrix"),
56
+setReplaceMethod("sigma2B", signature=signature(object="CNSet", value="ff_or_matrix"),
58 57
 	  function(object, value){
59 58
 		  linearParamElementReplace(object, "sig2B", value)
60 59
 })
61 60
 
62
-setReplaceMethod("tau2A", signature=signature(object="CNSet", value="matrix"),
61
+setReplaceMethod("tau2A", signature=signature(object="CNSet", value="ff_or_matrix"),
63 62
 	  function(object, value){
64 63
 		  linearParamElementReplace(object, "tau2A", value)
65 64
 })
66 65
 
67
-setReplaceMethod("tau2B", signature=signature(object="CNSet", value="matrix"),
66
+setReplaceMethod("tau2B", signature=signature(object="CNSet", value="ff_or_matrix"),
68 67
 	  function(object, value){
69 68
 		  linearParamElementReplace(object, "tau2B", value)
70 69
 })
71 70
 
72
-setReplaceMethod("corrAA", signature=signature(object="CNSet", value="matrix"),
71
+setReplaceMethod("corrAA", signature=signature(object="CNSet", value="ff_or_matrix"),
73 72
 	  function(object, value){
74 73
 		  linearParamElementReplace(object, "corrAA", value)
75 74
 })
76 75
 
77
-setReplaceMethod("corrAB", signature=signature(object="CNSet", value="matrix"),
76
+setReplaceMethod("corrAB", signature=signature(object="CNSet", value="ff_or_matrix"),
78 77
 	  function(object, value){
79 78
 		  linearParamElementReplace(object, "corrAB", value)
80 79
 })
81 80
 
82
-setReplaceMethod("corrBB", signature=signature(object="CNSet", value="matrix"),
81
+setReplaceMethod("corrBB", signature=signature(object="CNSet", value="ff_or_matrix"),
83 82
 	  function(object, value){
84 83
 		  linearParamElementReplace(object, "corrBB", value)
85 84
 })
86 85
 
86
+setReplaceMethod("flags", signature=signature(object="CNSet", value="ff_or_matrix"),
87
+		 function(object, value){
88
+			 linearParamElementReplace(object, "flags", value)
89
+})
90
+##setReplaceMethod("flags", signature=signature(object="CNSet", value="ff_matrix"),
91
+##		 function(object, value){
92
+##			 linearParamElementReplace(object, "flags", value)
93
+##})
94
+
95
+ setReplaceMethod("numberGenotype",
96
+		  signature=signature(object="CNSet", value="ff_or_matrix"),
97
+		 function(object, value, batchname){
98
+			 assayDataElementReplace(object, batchname, value)
99
+		 })
100
+
87 101
 
88 102
 ##setMethod("ellipse", "CNSet", function(x, copynumber, batch, ...){
89 103
 ##	ellipse.CNSet(x, copynumber, batch, ...)
90 104
deleted file mode 100644
... ...
@@ -1 +0,0 @@
1
-
... ...
@@ -107,7 +107,7 @@ batch variable.
107 107
 <<celfiles>>=
108 108
 celFiles <- list.celfiles(pathToCels, full.names=TRUE, pattern=".CEL")
109 109
 celFiles <- celFiles[substr(basename(celFiles), 13, 13) %in% c("C", "Y")]
110
-batch <- substr(basename(celFiles), 13, 13)
110
+batch <- as.factor(substr(basename(celFiles), 13, 13))
111 111
 @
112 112
 
113 113
 The preprocessing steps for copy number estimation includes quantile
... ...
@@ -141,14 +141,38 @@ file implies that we have already run the copy number estimation and,
141 141
 therefore, we do not need to preprocess and genotype the files a
142 142
 second time.
143 143
 
144
-<<genotype, eval=FALSE>>=
144
+<<genotype>>=
145
+tmpdir <- "~/tmp"
146
+stop()
147
+if(FALSE){
148
+if(file.exists(tmpdir)) unlink(tmpdir, recursive=TRUE)
149
+dir.create(tmpdir)
150
+##trace(genotype, browser)
151
+system.time(gtSet2 <- genotype(filenames=celFiles[1:75],
152
+			       cdfName=cdfName,
153
+			       batch=batch[1:75]))
154
+save(gtSet2, file=file.path(tmpdir, "gtSet2.rda"))
155
+library(ff)
156
+ldPath(tmpdir)
157
+system.time(gtSet <- genotype(filenames=celFiles[1:75],
158
+			  cdfName=cdfName,
159
+			  batch=batch[1:75]))
160
+save(gtSet, file=file.path(tmpdir, "gtSet.rda"))
161
+q("no")
162
+}
163
+load(file.path(tmpdir, "gtSet.rda"))
164
+load(file.path(tmpdir, "gtSet2.rda"))
165
+
166
+
167
+
168
+##using matrices
145 169
 if(!file.exists(file.path(outdir, "cnSet.assayData_matrix.rda"))){
146
-	gtSet.assayData_matrix <- checkExists("gtSet.assayData_matrix",
147
-					   .path=outdir,
148
-					   .FUN=genotype,
149
-					   filenames=celFiles[1:20],
150
-					   cdfName=cdfName,
151
-					   batch=batch[1:20])
170
+	gtSet <- checkExists("gtSet.assayData_matrix",
171
+			     .path=outdir,
172
+			     .FUN=genotype,
173
+			     filenames=celFiles[1:20],
174
+			     cdfName=cdfName,
175
+			     batch=batch[1:20])
152 176
 	class(calls(gtSet.assayData_matrix))
153 177
 }
154 178
 ##q("no")
... ...
@@ -31,6 +31,45 @@ data(sample.CNSetLM)
31 31
 cnSet <- as(sample.CNSetLM, "CNSet")
32 32
 all(isCurrent(cnSet)) ## is the cnSet object current?
33 33
 
34
+library(crlmm)
35
+library(ff)
36
+outdir <- "~/madman/Rpacks/crlmm/inst/extdata"
37
+ldPath(outdir)
38
+data(sample.CNSetLMff)
39
+cnSetff <- as(sample.CNSetLMff, "CNSet")
40
+all(isCurrent(cnSetff))
41
+
42
+## a bigger object with multiple batches
43
+\dontrun{
44
+outdir <- "/amber1/scratch/rscharpf/jss/hapmap2"
45
+load(file.path(outdir, "container.rda"))
46
+container <- object; rm(object)
47
+container2 <- as(container, "CNSet")
48
+all(isCurrent(container2))
49
+## test replacement methods, subset methods
50
+table(batch(container2))
51
+##generates warning ... would need open, close in the '[' method
52
+invisible(open(nuA(container2)))
53
+xx <- nu(container2, "A")[1:5, ]
54
+nuA(container2)[1:5, ] <- xx
55
+invisible(close(nuA(container2)))
56
+}
57
+
58
+
59
+
60
+
61
+nms <- names(lM(sample.CNSetLMff))
62
+##names must be tau2A, tau2B, ...
63
+nms2 <- sapply(nms, function(x) strsplit(x, "\\.1")[[1]][[1]], USE.NAMES=FALSE)
64
+names(sample.CNSetLMff@lM) <- nms2
65
+cnSetff <- as(sample.CNSetLMff, "CNSet")
66
+## try replacement method without warning
67
+nuA(cnSetff)[1:10, ] <- matrix(1:10, 10, 1)
68
+
69
+## try subset method without warning
70
+
71
+
72
+
34 73
 ## information on the features
35 74
 fvarLabels(cnSet)
36 75
 ##