Browse code

Added help files for copy number accessors. Fixed bugs in CA, CB, and total copy number methods.

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

Rob Scharp authored on 21/08/2010 02:49:14
Showing11 changed files

... ...
@@ -34,6 +34,7 @@ importMethodsFrom(oligoClasses, allele, calls, "calls<-", confs,
34 34
 importFrom(oligoClasses, chromosome2integer, celfileDate, list.celfiles,
35 35
            copyNumber, initializeBigMatrix, initializeBigVector)
36 36
 
37
+
37 38
 importFrom(graphics, abline, axis, layout, legend, mtext, par, plot,
38 39
            polygon, rect, segments, text, points, boxplot, lines)
39 40
 
... ...
@@ -57,12 +58,12 @@ importFrom(ff, ffdf, physical.ff, physical.ffdf)
57 58
 importClassesFrom(oligoClasses, ffdf)
58 59
 
59 60
 exportMethods(lines)
60
-exportMethods(CA, CB, totalCopyNumber)
61
-export(crlmm, 
62
-       crlmmIllumina, 
61
+exportMethods(CA, CB)
62
+export(crlmm,
63
+       crlmmIllumina,
63 64
        crlmmIllumina2,
64 65
        ellipseCenters,
65
-       readIdatFiles, 
66
+       readIdatFiles,
66 67
        readIdatFiles2,
67 68
        snprma,
68 69
        snprma2,
... ...
@@ -74,6 +75,7 @@ export(constructIlluminaCNSet)
74 75
 export(computeCN, fit.lm1, fit.lm2, fit.lm3, fit.lm4, construct,
75 76
        dqrlsWrapper, fit.wls)
76 77
 export(computeCopynumber, ACN)
78
+export(totalCopyNumber)
77 79
 
78 80
 ## For debugging
79 81
 ## exportPattern("^[^\\.]")
... ...
@@ -2,5 +2,4 @@ setOldClass("ellipse")
2 2
 setOldClass("ff_matrix")
3 3
 setOldClass("ffdf")
4 4
 setClassUnion("ff_or_matrix", c("ff_matrix", "matrix", "ffdf"))
5
-setClassUnion("integerOrMissing", c("integer", "missing", "numeric"))
6 5
 
... ...
@@ -7,11 +7,9 @@ setGeneric("computeCopynumber", function(object, ...) standardGeneric("computeCo
7 7
 setGeneric("snpIndex", function(object) standardGeneric("snpIndex"))
8 8
 setGeneric("snpNames", function(object) standardGeneric("snpNames"))
9 9
 
10
-setGeneric("totalCopyNumber", function(object,...) standardGeneric("totalCopyNumber"))
11
-
12
-setGeneric("CA", function(object, i, j, ...) standardGeneric("CA"))
13
-setGeneric("CB", function(object, i, j, ...) standardGeneric("CB"))
14
-setGeneric("totalCopyNumber", function(object, i, j, ...) standardGeneric("totalCopyNumber"))
10
+setGeneric("CA", function(object, ...) standardGeneric("CA"))
11
+setGeneric("CB", function(object, ...) standardGeneric("CB"))
12
+##setGeneric("totalCopyNumber", function(object, ...) standardGeneric("totalCopyNumber"))
15 13
 
16 14
 
17 15
 ## The generics below are for internal use with copy number methods
... ...
@@ -58,17 +58,6 @@ construct <- function(filenames, cdfName, copynumber=FALSE,
58 58
 	}
59 59
 	if(missing(sns) & missing(filenames)) stop("one of filenames or samplenames (sns) must be provided")
60 60
 	if(verbose) message("Initializing container for assay data elements alleleA, alleleB, call, callProbability")
61
-	if(!missing(filenames)){
62
-		if(missing(sns)) sns <- basename(filenames)
63
-		protocolData <- getProtocolData.Affy(filenames)
64
-		protocolData$batch <- as.numeric(as.factor(protocolData$ScanDate))
65
-	} else{
66
-		protocolData <- new("AnnotatedDataFrame",
67
-				    data=data.frame(batch=batch),
68
-				    varMetadata=data.frame(labelDescription="batch",
69
-				    row.names="batch"))
70
-	}
71
-	rownames(pData(protocolData)) <- sns
72 61
 	featureData <- getFeatureData.Affy(cdfName, copynumber=copynumber)
73 62
 	if(!missing(fns)){
74 63
 		index <- match(fns, featureNames(featureData))
... ...
@@ -76,25 +65,27 @@ construct <- function(filenames, cdfName, copynumber=FALSE,
76 65
 		featureData <- featureData[index, ]
77 66
 	}
78 67
 	nr <- nrow(featureData); nc <- length(sns)
79
-	adObjects <- list(alleleA=initializeBigMatrix(name="A", nr, nc),
80
-			  alleleB=initializeBigMatrix(name="B", nr, nc),
81
-			  call=initializeBigMatrix(name="call", nr, nc),
82
-			  callProbability=initializeBigMatrix(name="callPr", nr,nc))
68
+	cnSet <- new("CNSet",
69
+		     alleleA=initializeBigMatrix(name="A", nr, nc),
70
+		     alleleB=initializeBigMatrix(name="B", nr, nc),
71
+		     call=initializeBigMatrix(name="call", nr, nc),
72
+		     callProbability=initializeBigMatrix(name="callPr", nr,nc),
73
+		     annotation=cdfName,
74
+		     batch=batch)
75
+	sampleNames(cnSet) <- sns
76
+	if(!missing(filenames)){
77
+		if(missing(sns)) sns <- basename(filenames)
78
+		protocolData <- getProtocolData.Affy(filenames)
79
+	} else{
80
+		protocolData <- annotatedDataFrameFrom(A(cnSet), byrow=FALSE)
81
+	}
82
+	rownames(pData(protocolData)) <- sns
83
+	protocolData(callSet) <- protocolData
84
+	featureData(cnSet) <- featureData
83 85
 	pd <- data.frame(matrix(NA, nc, 3), row.names=sns)
84 86
 	colnames(pd)=c("SKW", "SNR", "gender")
85
-	phenoData <- new("AnnotatedDataFrame", data=pd)
86
-	adObjects <- lapply(adObjects, function(x,sns) {colnames(x) <- sns; return(x)}, sns=sns)
87
-	callSet <- new("CNSet",
88
-		       alleleA=adObjects[["alleleA"]],
89
-		       alleleB=adObjects[["alleleB"]],
90
-		       call=adObjects[["call"]],
91
-		       callProbability=adObjects[["callProbability"]],
92
-		       protocolData=protocolData,
93
-		       phenoData=phenoData,
94
-		       featureData=featureData,
95
-		       annotation=cdfName)
96
-	lM(callSet) <- initializeParamObject(list(featureNames(callSet), unique(protocolData(callSet)$batch)))
97
-	return(callSet)
87
+	phenoData(cnSet) <- new("AnnotatedDataFrame", data=pd)
88
+	return(cnSet)
98 89
 }
99 90
 
100 91
 ##genotype <- function(filenames,
... ...
@@ -391,7 +382,7 @@ fit.wls <- function(allele, Ystar, W, Ns, autosome=TRUE){
391 382
 		if(allele == "B") X <- cbind(1, c(0, 1, 0, 1, 2), c(1, 0, 2, 1, 0))
392 383
 		betahat <- matrix(NA, 3, nrow(Ystar))
393 384
 	}
394
-	if(NOHET) X <- X[-2, ] ##more than 1 X chromosome, but all homozygous		
385
+	if(NOHET) X <- X[-2, ] ##more than 1 X chromosome, but all homozygous
395 386
 	##How to quickly generate Xstar, Xstar = diag(W) %*% X
396 387
 	##Xstar <- apply(W, 1, generateX, X)
397 388
 	ww <- rep(1, ncol(Ystar))
... ...
@@ -403,7 +394,7 @@ fit.wls <- function(allele, Ystar, W, Ns, autosome=TRUE){
403 394
 }
404 395
 
405 396
 ##nuphiAlleleX <- function(allele, Ystar, W, Ns, chrX=FALSE){
406
-##	complete <- rowSums(is.na(W)) == 0 
397
+##	complete <- rowSums(is.na(W)) == 0
407 398
 ##	if(any(!is.finite(W))){## | any(!is.finite(V))){
408 399
 ##		i <- which(rowSums(!is.finite(W)) > 0)
409 400
 ##		stop("Possible zeros in the within-genotype estimates of the spread (vA, vB). ")
... ...
@@ -411,14 +402,14 @@ fit.wls <- function(allele, Ystar, W, Ns, autosome=TRUE){
411 402
 ##	##NOHET <- mean(Ns[, 2], na.rm=TRUE) < 0.05
412 403
 ##	if(missing(allele)) stop("must specify allele")
413 404
 ##	if(allele == "A") X <- cbind(1, c(1, 0, 2, 1, 0), c(0, 1, 0, 1, 2))
414
-##	if(allele == "B") X <- cbind(1, c(0, 1, 0, 1, 2), c(1, 0, 2, 1, 0))	
415
-##	##if(NOHET) X <- X[-2, ] ##more than 1 X chromosome, but all homozygous		
405
+##	if(allele == "B") X <- cbind(1, c(0, 1, 0, 1, 2), c(1, 0, 2, 1, 0))
406
+##	##if(NOHET) X <- X[-2, ] ##more than 1 X chromosome, but all homozygous
416 407
 ##	##How to quickly generate Xstar, Xstar = diag(W) %*% X
417 408
 ##	Xstar <- apply(W, 1, generateX, X)
418 409
 ##	IXTX <- apply(Xstar, 2, generateIXTX, nrow=nrow(X))
419 410
 ##
420 411
 ##	betahat <- matrix(NA, 3, nrow(Ystar))
421
-##	ses <- matrix(NA, 3, nrow(Ystar))			
412
+##	ses <- matrix(NA, 3, nrow(Ystar))
422 413
 ##	##betahat <- matrix(NA, 2, nrow(Ystar))
423 414
 ##	##ses <- matrix(NA, 2, nrow(Ystar))
424 415
 ##	for(i in 1:nrow(Ystar)){
... ...
@@ -442,8 +433,8 @@ fit.wls <- function(allele, Ystar, W, Ns, autosome=TRUE){
442 433
 ##	if(any(!is.finite(W))){## | any(!is.finite(V))){
443 434
 ##		i <- which(rowSums(!is.finite(W)) > 0)
444 435
 ##		stop("Possible zeros in the within-genotype estimates of the spread (vA, vB). ")
445
-##	}	
446
-##	Ns <- tmp.objects[["Ns"]]	
436
+##	}
437
+##	Ns <- tmp.objects[["Ns"]]
447 438
 ##	Ns <- Ns[I, ]
448 439
 ##	CHR <- unique(chromosome(object))
449 440
 ##	batch <- unique(batch(object))
... ...
@@ -453,7 +444,7 @@ fit.wls <- function(allele, Ystar, W, Ns, autosome=TRUE){
453 444
 ##	phiB <- getParam(object, "phiB", batch)
454 445
 ##	if(CHR==23){
455 446
 ##		phiAX <- getParam(object, "phiAX", batch)
456
-##		phiBX <- getParam(object, "phiBX", batch)		
447
+##		phiBX <- getParam(object, "phiBX", batch)
457 448
 ##	}
458 449
 ##	NOHET <- mean(Ns[, "AB"], na.rm=TRUE) < 0.05
459 450
 ##	if(missing(allele)) stop("must specify allele")
... ...
@@ -483,14 +474,14 @@ fit.wls <- function(allele, Ystar, W, Ns, autosome=TRUE){
483 474
 ##		object <- pr(object, "phiA", batch, phiA)
484 475
 ##		if(CHR == 23){
485 476
 ##			phiAX[complete] <- betahat[3, ]
486
-##			object <- pr(object, "phiAX", batch, phiAX)			
477
+##			object <- pr(object, "phiAX", batch, phiAX)
487 478
 ##		}
488 479
 ##	}
489 480
 ##	if(allele=="B"){
490 481
 ##		nuB[complete] <- betahat[1, ]
491 482
 ##		phiB[complete] <- betahat[2, ]
492 483
 ##		object <- pr(object, "nuB", batch, nuB)
493
-##		object <- pr(object, "phiB", batch, phiB)		
484
+##		object <- pr(object, "phiB", batch, phiB)
494 485
 ##		if(CHR == 23){
495 486
 ##			phiBX[complete] <- betahat[3, ]
496 487
 ##			object <- pr(object, "phiBX", batch, phiBX)
... ...
@@ -507,7 +498,7 @@ predictGender <- function(res, cdfName="genomewidesnp6", SNRMin=5){
507 498
         loader("23index.rda", .crlmmPkgEnv, pkgname)
508 499
 	index <- getVarInEnv("index")
509 500
 	##load(file.path(path, "23index.rda"))
510
-	XIndex <- index[[1]]	
501
+	XIndex <- index[[1]]
511 502
 	if(class(res) != "ABset"){
512 503
 		A <- res$A
513 504
 		B <- res$B
... ...
@@ -629,7 +620,7 @@ predictGender <- function(res, cdfName="genomewidesnp6", SNRMin=5){
629 620
 ##					column <- match(batchName, colnames(lM(object)[[k]]))
630 621
 ##					lM(object)[[k]][row.index, column] <- fData(tmp)[, k]
631 622
 ##				}
632
-##			}			
623
+##			}
633 624
 ##			rm(tmp); ##gc()
634 625
 ##			ii <- ii+1
635 626
 ##		}
... ...
@@ -666,7 +657,7 @@ crlmmCopynumberLD <- function(object,
666 657
 	##autosomeIndex.snps <- (1:nrow(object))[chromosome(object) < 23 & isSnp(object) & !is.na(chromosome(object))]
667 658
 	autosomeIndex.snps <- which(chromosome(object) < 23 & isSnp(object) & !is.na(chromosome(object)))
668 659
 	autosomeIndex.nps <- which(chromosome(object) < 23 & !isSnp(object) & !is.na(chromosome(object)))
669
-	##autosomeIndex.nps <- (1:nrow(object))[chromosome(object) < 23 & !isSnp(object) & !is.na(chromosome(object))]	
660
+	##autosomeIndex.nps <- (1:nrow(object))[chromosome(object) < 23 & !isSnp(object) & !is.na(chromosome(object))]
670 661
 
671 662
 	## Do chromosome X in batches
672 663
 	Ns <- initializeBigMatrix("Ns", nrow(object), 5)
... ...
@@ -682,7 +673,7 @@ crlmmCopynumberLD <- function(object,
682 673
 		save(snpflags, file=file.path(ldPath(), "snpflags.rda"))
683 674
 	} else{
684 675
 		load(file.path(ldPath(), "snpflags.rda"))
685
-	} 
676
+	}
686 677
 	if(verbose) message("Estimating allele-specific copy number at autosomal SNPs")
687 678
 	snpBatches <- splitIndicesByLength(autosomeIndex.snps, ocProbesets())
688 679
 	ocLapply(seq(along=snpBatches),
... ...
@@ -795,7 +786,7 @@ fit.lm1 <- function(idxBatch,
795 786
 	snps <- snpBatches[[idxBatch]]
796 787
 	batches <- split(seq(along=batch(object)), batch(object))
797 788
 	batches <- batches[sapply(batches, length) >= MIN.SAMPLES]
798
-	
789
+
799 790
 	open(object)
800 791
 	open(snpflags)
801 792
 	open(normal)
... ...
@@ -904,14 +895,14 @@ fit.lm1 <- function(idxBatch,
904 895
 		res <- fit.wls(allele="A", Ystar=YA, W=wA, Ns=Ns)
905 896
 ##		} else{
906 897
 ##			if(zzzz==1) message("currently, only weighted least squares (wls) is available... fitting wls")
907
-##			res <- fit.wls(allele="A", Ystar=YA, W=wA, Ns=Ns)			
898
+##			res <- fit.wls(allele="A", Ystar=YA, W=wA, Ns=Ns)
908 899
 ##		}
909 900
 		nuA[, J] <- res[[1]]
910 901
 		phiA[, J] <- res[[2]]
911 902
 		res <- fit.wls(allele="B", Ystar=YB, W=wB, Ns=Ns)
912 903
 ##		} else {
913 904
 ##			if(zzzz==1) message("currently, only weighted least squares (wls) is available... fitting wls")
914
-##			res <- fit.wls(allele="B", Ystar=YB, W=wB, Ns=Ns)			
905
+##			res <- fit.wls(allele="B", Ystar=YB, W=wB, Ns=Ns)
915 906
 ##		}
916 907
 		##nuB[, J] <- res[[1]]
917 908
 		nuB[, J] <- res[1, ]
... ...
@@ -939,7 +930,7 @@ fit.lm1 <- function(idxBatch,
939 930
 
940 931
 	snpflags[snps, ] <- flags
941 932
 	lapply(lM(object), open)
942
-	
933
+
943 934
 	tmp <- physical(lM(object))$tau2A
944 935
 	tmp[snps, ] <- tau2A
945 936
 	lM(object)$tau2A <- tmp
... ...
@@ -976,7 +967,7 @@ fit.lm1 <- function(idxBatch,
976 967
 	tmp <- physical(lM(object))$corrBB
977 968
 	tmp[snps, ] <- corrBB
978 969
 	lM(object)$corrBB <- tmp
979
-	
970
+
980 971
 	lapply(assayData(object), close)
981 972
 	lapply(lM(object), close)
982 973
 	TRUE
... ...
@@ -1004,14 +995,14 @@ fit.lm2 <- function(idxBatch,
1004 995
 	snps <- snpBatches[[idxBatch]]
1005 996
 	batches <- split(seq(along=batch(object)), batch(object))
1006 997
 	batches <- batches[sapply(batches, length) >= MIN.SAMPLES]
1007
-	
998
+
1008 999
 	open(object)
1009 1000
 	open(snpflags)
1010 1001
 	open(normal)
1011 1002
 
1012 1003
 ##	cA <- matrix(NA, length(snps), ncol(object))
1013 1004
 	ii <- isSnp(object) & chromosome(object) < 23 & !is.na(chromosome(object))
1014
-	flags <- as.matrix(snpflags[,])  
1005
+	flags <- as.matrix(snpflags[,])
1015 1006
 	noflags <- rowSums(flags, na.rm=TRUE) == 0  ##NA's for unevaluated batches
1016 1007
 	## We do not want to write to discuss for each batch.  More efficient to
1017 1008
 	## write to disk after estimating these parameters for all batches.
... ...
@@ -1115,7 +1106,7 @@ fit.lm3 <- function(idxBatch,
1115 1106
 		snps <- snpBatches[[idxBatch]]
1116 1107
 	batches <- split(seq(along=batch(object)), batch(object))
1117 1108
 	batches <- batches[sapply(batches, length) >= MIN.SAMPLES]
1118
-	
1109
+
1119 1110
 	open(snpflags)
1120 1111
 	open(normal)
1121 1112
 	open(object)
... ...
@@ -1133,7 +1124,7 @@ fit.lm3 <- function(idxBatch,
1133 1124
 	AA <- as.matrix(A(object)[snps, ])
1134 1125
 	BB <- as.matrix(B(object)[snps, ])
1135 1126
 	for(k in batches){
1136
-		##if(verbose) message("SNP batch ", ii, " of ", length(batches)) 
1127
+		##if(verbose) message("SNP batch ", ii, " of ", length(batches))
1137 1128
 		## within-genotype moments
1138 1129
 		gender <- object$gender[k]
1139 1130
 		G <- GG[, k]
... ...
@@ -1154,7 +1145,7 @@ fit.lm3 <- function(idxBatch,
1154 1145
 		vA.F <- applyByGenotype(A[, gender==2], rowMAD, G[, gender==2])
1155 1146
 		vB.F <- applyByGenotype(B[, gender==2], rowMAD, G[, gender==2])
1156 1147
 		vA.M <- applyByGenotype(A[, gender==1], rowMAD, G[, gender==1])
1157
-		vB.M <- applyByGenotype(B[, gender==1], rowMAD, G[, gender==1])		
1148
+		vB.M <- applyByGenotype(B[, gender==1], rowMAD, G[, gender==1])
1158 1149
 		vA.F <- shrink(vA.F, Ns.F, DF.PRIOR)
1159 1150
 		vA.M <- shrink(vA.M, Ns.M, DF.PRIOR)
1160 1151
 		vB.F <- shrink(vB.F, Ns.F, DF.PRIOR)
... ...
@@ -1197,7 +1188,7 @@ fit.lm3 <- function(idxBatch,
1197 1188
 		notMissing <- !(is.na(muA.M[, 1]) | is.na(muA.M[, 3]) | is.na(muB.M[, 1]) | is.na(muB.M[, 3]))
1198 1189
 		complete <- list()
1199 1190
 		complete[[1]] <- which(correct.orderA & correct.orderB & nobsA & notMissing) ##be selective here
1200
-		complete[[2]] <- which(correct.orderA & correct.orderB & nobsB & notMissing) ##be selective here	
1191
+		complete[[2]] <- which(correct.orderA & correct.orderB & nobsB & notMissing) ##be selective here
1201 1192
 		size <- min(5000, length(complete[[1]]))
1202 1193
 		if(size > 5000) complete <- lapply(complete, function(x) sample(x, size))
1203 1194
 		##
... ...
@@ -1273,7 +1264,7 @@ fit.lm3 <- function(idxBatch,
1273 1264
 ##	cB[cB < 0.05] <- 0.05
1274 1265
 ##	cA[cA > 5] <-  5
1275 1266
 ##	cB[cB > 5] <- 5
1276
-	
1267
+
1277 1268
 	##--------------------------------------------------
1278 1269
 	##RS: need to fix.  why are there NA's by coercion
1279 1270
 ##	cA <- matrix(as.integer(cA*100), nrow(cA), ncol(cA))
... ...
@@ -1317,7 +1308,7 @@ fit.lm3 <- function(idxBatch,
1317 1308
 	tmp <- physical(lM(object))$phiPrimeB
1318 1309
 	tmp[snps, ] <- phiB2
1319 1310
 	lM(object)$phiPrimeB <- tmp
1320
-	
1311
+
1321 1312
 	tmp <- physical(lM(object))$corrAB
1322 1313
 	tmp[snps, ] <- corrAB
1323 1314
 	lM(object)$corrAB <- tmp
... ...
@@ -1326,7 +1317,7 @@ fit.lm3 <- function(idxBatch,
1326 1317
 	lM(object)$corrAA <- tmp
1327 1318
 	tmp <- physical(lM(object))$corrBB
1328 1319
 	tmp[snps, ] <- corrBB
1329
-	lM(object)$corrBB <- tmp	
1320
+	lM(object)$corrBB <- tmp
1330 1321
 	lapply(assayData(object), close)
1331 1322
 	lapply(lM(object), close)
1332 1323
 	TRUE
... ...
@@ -1350,7 +1341,7 @@ fit.lm4 <- function(idxBatch,
1350 1341
 		    MIN.PHI,
1351 1342
 		    verbose, ...){
1352 1343
 	physical <- get("physical")
1353
-	if(verbose) message("Probe stratum ", idxBatch, " of ", length(snpBatches))	
1344
+	if(verbose) message("Probe stratum ", idxBatch, " of ", length(snpBatches))
1354 1345
 	open(object)
1355 1346
 	open(normal)
1356 1347
 	open(snpflags)
... ...
@@ -1372,7 +1363,7 @@ fit.lm4 <- function(idxBatch,
1372 1363
 	i2 <- rowSums(nuIB < 20, na.rm=TRUE) == 0
1373 1364
 	i3 <- rowSums(phiIA < 20, na.rm=TRUE) == 0
1374 1365
 	i4 <- rowSums(phiIB < 20, na.rm=TRUE) == 0
1375
-	
1366
+
1376 1367
 	snp.index <- which(i1 & i2 & i3 & i4 & noflags)
1377 1368
 	if(length(snp.index) == 0){
1378 1369
 		warning("No snps meet the following criteria: (1) nu and phi > 20 and (2) at least MIN.OBS in each genotype cluster. CN not estimated for nonpolymorphic loci on X")
... ...
@@ -1399,7 +1390,7 @@ fit.lm4 <- function(idxBatch,
1399 1390
 	##if(missing(which.batches)) which.batches <- seq(along=batches)
1400 1391
 	##batches <- batches[which.batches]
1401 1392
 	for(k in batches){
1402
-		##if(verbose) message("SNP batch ", ii, " of ", length(batches)) 
1393
+		##if(verbose) message("SNP batch ", ii, " of ", length(batches))
1403 1394
 		G <- GG[, k]
1404 1395
 		xx <- CP[, k]
1405 1396
 		highConf <- (1-exp(-xx/1000)) > GT.CONF.THR
... ...
@@ -1408,7 +1399,7 @@ fit.lm4 <- function(idxBatch,
1408 1399
 		AA <- A.snp[, k]
1409 1400
 		BB <- B.snp[, k]
1410 1401
 
1411
-				  
1402
+
1412 1403
 		##index <- GT.B <- GT.A <- vector("list", 3)
1413 1404
 		##names(index) <- names(GT.B) <- names(GT.A) <- c("AA", "AB", "BB")
1414 1405
 		Ns <- applyByGenotype(matrix(1, nrow(G), ncol(G)), rowSums, G)
... ...
@@ -1418,7 +1409,7 @@ fit.lm4 <- function(idxBatch,
1418 1409
 		muB <- muB[, 3]
1419 1410
 		X <- cbind(1, log2(c(muA, muB)))
1420 1411
 		J <- match(unique(batch(object)[k]), unique(batch(object)))
1421
-		
1412
+
1422 1413
 		Y <- log2(c(phiA.snp[, J], phiB.snp[, J]))
1423 1414
 
1424 1415
 		##--------------------------------------------------
... ...
@@ -1428,7 +1419,7 @@ fit.lm4 <- function(idxBatch,
1428 1419
 		X <- X[!remove, ]
1429 1420
 		##--------------------------------------------------
1430 1421
 		betahat <- solve(crossprod(X), crossprod(X, Y))
1431
-		
1422
+
1432 1423
 
1433 1424
 		##nonpolymorphic
1434 1425
 		A <- AA.np[, k]
... ...
@@ -1522,7 +1513,7 @@ cnrma2 <- function(A, filenames, row.names, verbose=TRUE, seed=1, cdfName, sns){
1522 1513
 		 row.names=row.names,
1523 1514
 		 filenames=filenames,
1524 1515
 		 A=A,
1525
-		 SKW=SKW, 
1516
+		 SKW=SKW,
1526 1517
 		 seed=seed,
1527 1518
 		 pkgname=pkgname,
1528 1519
 		 cdfName=cdfName,
... ...
@@ -1656,15 +1647,15 @@ processCEL2 <- function(i, filenames, row.names, A, SKW, seed, cdfName, pkgname)
1656 1647
 ##			    paste("sig2A", uplate, sep="_"),
1657 1648
 ##			    paste("sig2B", uplate, sep="_"),
1658 1649
 ##			    paste("nuA", uplate, sep="_"),
1659
-##			    paste("nuA.se", uplate, sep="_"),			    
1650
+##			    paste("nuA.se", uplate, sep="_"),
1660 1651
 ##			    paste("nuB", uplate, sep="_"),
1661
-##			    paste("nuB.se", uplate, sep="_"),			    			    
1652
+##			    paste("nuB.se", uplate, sep="_"),
1662 1653
 ##			    paste("phiA", uplate, sep="_"),
1663
-##			    paste("phiA.se", uplate, sep="_"),			    
1654
+##			    paste("phiA.se", uplate, sep="_"),
1664 1655
 ##			    paste("phiB", uplate, sep="_"),
1665
-##			    paste("phiB.se", uplate, sep="_"),			    
1656
+##			    paste("phiB.se", uplate, sep="_"),
1666 1657
 ##			    paste("phiAX", uplate, sep="_"),
1667
-##			    paste("phiBX", uplate, sep="_"),			    
1658
+##			    paste("phiBX", uplate, sep="_"),
1668 1659
 ##			    paste("corr", uplate, sep="_"),
1669 1660
 ##			    paste("corrA.BB", uplate, sep="_"),
1670 1661
 ##			    paste("corrB.AA", uplate, sep="_"))
... ...
@@ -1703,7 +1694,7 @@ processCEL2 <- function(i, filenames, row.names, A, SKW, seed, cdfName, pkgname)
1703 1694
 ##	nuA <- getParam(object, "nuA", batch)
1704 1695
 ##	nuB <- getParam(object, "nuB", batch)
1705 1696
 ##	phiA <- getParam(object, "phiA", batch)
1706
-##	phiB <- getParam(object, "phiB", batch)	
1697
+##	phiB <- getParam(object, "phiB", batch)
1707 1698
 ##	sns <- sampleNames(object)
1708 1699
 ##	muA <- tmp.objects[["muA"]]
1709 1700
 ##	muB <- tmp.objects[["muB"]]
... ...
@@ -1733,13 +1724,13 @@ processCEL2 <- function(i, filenames, row.names, A, SKW, seed, cdfName, pkgname)
1733 1724
 ##		ix <- 1:nrow(X)
1734 1725
 ##	}
1735 1726
 ##	betahat <- solve(crossprod(X[ix, ]), crossprod(X[ix, ], Y[ix]))
1736
-##	normal <- tmp.objects[["normal"]][!isSnp(object), ]	
1727
+##	normal <- tmp.objects[["normal"]][!isSnp(object), ]
1737 1728
 ##	if(CHR == 23){
1738 1729
 ##		##normalNP <- envir[["normalNP"]]
1739 1730
 ##		##normalNP <- normalNP[, plate==uplate[p]]
1740 1731
 ##		##nuT <- envir[["nuT"]]
1741 1732
 ##		##phiT <- envir[["phiT"]]
1742
-##		
1733
+##
1743 1734
 ##		##cnvs <- envir[["cnvs"]]
1744 1735
 ##                ##loader("cnProbes.rda", pkgname=pkgname, envir=.crlmmPkgEnv)
1745 1736
 ##                ##cnProbes <- get("cnProbes", envir=.crlmmPkgEnv)
... ...
@@ -1767,7 +1758,7 @@ processCEL2 <- function(i, filenames, row.names, A, SKW, seed, cdfName, pkgname)
1767 1758
 ##		mus <- log2(cbind(mu1, mu2))
1768 1759
 ##		X.men <- cbind(1, mus[, 1])
1769 1760
 ##		X.fem <- cbind(1, mus[, 2])
1770
-##		
1761
+##
1771 1762
 ##		Yhat1 <- as.numeric(X.men %*% betahat)
1772 1763
 ##		Yhat2 <- as.numeric(X.fem %*% betahat)
1773 1764
 ##		phi1 <- 2^(Yhat1)
... ...
@@ -1796,17 +1787,17 @@ processCEL2 <- function(i, filenames, row.names, A, SKW, seed, cdfName, pkgname)
1796 1787
 ##
1797 1788
 ##		nuA[!isSnp(object)] <- nu2
1798 1789
 ##		phiA[!isSnp(object)] <- phi2
1799
-##		
1790
+##
1800 1791
 ##		THR.NU.PHI <- cnOptions$THR.NU.PHI
1801 1792
 ##		if(THR.NU.PHI){
1802 1793
 ##			verbose <- cnOptions$verbose
1803 1794
 ##			##Assign values to object
1804 1795
 ##			object <- pr(object, "nuA", batch, nuA)
1805
-##			object <- pr(object, "phiA", batch, phiA)			
1796
+##			object <- pr(object, "phiA", batch, phiA)
1806 1797
 ##			##if(verbose) message("Thresholding nu and phi")
1807 1798
 ##			object <- thresholdModelParams(object, cnOptions)
1808 1799
 ##		} else {
1809
-##			object <- pr(object, "nuA", batch, nuA)		
1800
+##			object <- pr(object, "nuA", batch, nuA)
1810 1801
 ##			object <- pr(object, "phiA", batch, phiA)
1811 1802
 ##		}
1812 1803
 ##	} else {
... ...
@@ -1823,7 +1814,7 @@ processCEL2 <- function(i, filenames, row.names, A, SKW, seed, cdfName, pkgname)
1823 1814
 ##			verbose <- cnOptions$verbose
1824 1815
 ##			##Assign values to object
1825 1816
 ##			object <- pr(object, "nuA", batch, nuA)
1826
-##			object <- pr(object, "phiA", batch, phiA)			
1817
+##			object <- pr(object, "phiA", batch, phiA)
1827 1818
 ##			##if(verbose) message("Thresholding nu and phi")
1828 1819
 ##			object <- thresholdModelParams(object, cnOptions)
1829 1820
 ##			##reassign values (now thresholded at MIN.NU and MIN.PHI
... ...
@@ -1851,7 +1842,7 @@ processCEL2 <- function(i, filenames, row.names, A, SKW, seed, cdfName, pkgname)
1851 1842
 ##	vA <- tmp.objects[["vA"]]
1852 1843
 ##	vB <- tmp.objects[["vB"]]
1853 1844
 ##	Ns <- tmp.objects[["Ns"]]
1854
-##	G <- snpCall(object) 
1845
+##	G <- snpCall(object)
1855 1846
 ##	GT.CONF.THR <- cnOptions$GT.CONF.THR
1856 1847
 ##	CHR <- unique(chromosome(object))
1857 1848
 ##	A <- A(object)
... ...
@@ -1878,11 +1869,11 @@ processCEL2 <- function(i, filenames, row.names, A, SKW, seed, cdfName, pkgname)
1878 1869
 ##	for(j in 1:3){
1879 1870
 ##		GT <- G==j & highConf & IX & snpIndicator
1880 1871
 ##		GT <- GT * normal
1881
-##		Ns[, j+2] <- rowSums(GT, na.rm=TRUE)				
1872
+##		Ns[, j+2] <- rowSums(GT, na.rm=TRUE)
1882 1873
 ##		GT[GT == FALSE] <- NA
1883 1874
 ##		GT.A[[j]] <- GT*A
1884 1875
 ##		GT.B[[j]] <- GT*B
1885
-##		index[[j]] <- which(Ns[, j+2] > 0 & isSnp(object)) ##RS: added		
1876
+##		index[[j]] <- which(Ns[, j+2] > 0 & isSnp(object)) ##RS: added
1886 1877
 ##		muA[, j+2] <- rowMedians(GT.A[[j]], na.rm=TRUE)
1887 1878
 ##		muB[, j+2] <- rowMedians(GT.B[[j]], na.rm=TRUE)
1888 1879
 ##		vA[, j+2] <- rowMAD(GT.A[[j]], na.rm=TRUE)
... ...
@@ -1904,14 +1895,14 @@ processCEL2 <- function(i, filenames, row.names, A, SKW, seed, cdfName, pkgname)
1904 1895
 ##	if(CHR == 23){
1905 1896
 ##		k <- 1
1906 1897
 ##		for(j in c(1,3)){
1907
-##			GT <- G==j & highConf & !IX 
1898
+##			GT <- G==j & highConf & !IX
1908 1899
 ##			Ns[, k] <- rowSums(GT)
1909 1900
 ##			GT[GT == FALSE] <- NA
1910 1901
 ##			muA[, k] <- rowMedians(GT*A, na.rm=TRUE)
1911 1902
 ##			muB[, k] <- rowMedians(GT*B, na.rm=TRUE)
1912 1903
 ##			vA[, k] <- rowMAD(GT*A, na.rm=TRUE)
1913 1904
 ##			vB[, k] <- rowMAD(GT*B, na.rm=TRUE)
1914
-##			
1905
+##
1915 1906
 ##			DF <- Ns[, k]-1
1916 1907
 ##			DF[DF < 1] <- 1
1917 1908
 ##			v0A <- median(vA[, k], na.rm=TRUE)
... ...
@@ -1919,7 +1910,7 @@ processCEL2 <- function(i, filenames, row.names, A, SKW, seed, cdfName, pkgname)
1919 1910
 ##			vA[, k] <- (vA[, k]*DF + v0A*DF.PRIOR)/(DF.PRIOR+DF)
1920 1911
 ##			vA[is.na(vA[, k]), k] <- v0A
1921 1912
 ##			vB[, k] <- (vB[, k]*DF + v0B*DF.PRIOR)/(DF.PRIOR+DF)
1922
-##			vB[is.na(vB[, k]), k] <- v0B			
1913
+##			vB[is.na(vB[, k]), k] <- v0B
1923 1914
 ##			k <- k+1
1924 1915
 ##		}
1925 1916
 ##	}
... ...
@@ -1991,7 +1982,7 @@ imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
1991 1982
 ##	index.BB <- which(Ns[, "BB"] >= MIN.OBS)
1992 1983
 ##	correct.orderA <- muA[, "AA"] > muA[, "BB"]
1993 1984
 ##	correct.orderB <- muB[, "BB"] > muB[, "AA"]
1994
-##	##For chr X, this will ignore the males 
1985
+##	##For chr X, this will ignore the males
1995 1986
 ##	nobs <- rowSums(Ns[, 3:5] >= MIN.OBS, na.rm=TRUE) == 3
1996 1987
 ##	index.complete <- which(correct.orderA & correct.orderB & nobs) ##be selective here
1997 1988
 ##	size <- min(5000, length(index.complete))
... ...
@@ -2020,7 +2011,7 @@ imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
2020 2011
 ##	notMissing <- !(is.na(muA[, "A"]) | is.na(muA[, "B"]) | is.na(muB[, "A"]) | is.na(muB[, "B"]))
2021 2012
 ##	complete <- list()
2022 2013
 ##	complete[[1]] <- which(correct.orderA & correct.orderB & nobsA & notMissing) ##be selective here
2023
-##	complete[[2]] <- which(correct.orderA & correct.orderB & nobsB & notMissing) ##be selective here	
2014
+##	complete[[2]] <- which(correct.orderA & correct.orderB & nobsB & notMissing) ##be selective here
2024 2015
 ##	size <- min(5000, length(complete[[1]]))
2025 2016
 ##	if(size > 5000) complete <- lapply(complete, function(x) sample(x, size))
2026 2017
 ##	if(CHR == 23){
... ...
@@ -2070,7 +2061,7 @@ imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
2070 2061
 ##		muB[index[[j]], -c(1, 2, k+2)] <- mus[, 3:4]
2071 2062
 ##	}
2072 2063
 ##	negA <- rowSums(muA < 0) > 0
2073
-##	negB <- rowSums(muB < 0) > 0	
2064
+##	negB <- rowSums(muB < 0) > 0
2074 2065
 ##	snpflags <- snpflags| negA | negB | rowSums(is.na(muA[, 3:5]), na.rm=TRUE) > 0
2075 2066
 ##	tmp.objects$snpflags <- snpflags
2076 2067
 ##	tmp.objects[["muA"]] <- muA
... ...
@@ -2094,7 +2085,7 @@ imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
2094 2085
 ##	AA.A <- GT.A[[1]]
2095 2086
 ##	AB.A <- GT.A[[2]]
2096 2087
 ##	BB.A <- GT.A[[3]]
2097
-##	
2088
+##
2098 2089
 ##	AA.B <- GT.B[[1]]
2099 2090
 ##	AB.B <- GT.B[[2]]
2100 2091
 ##	BB.B <- GT.B[[3]]
... ...
@@ -2102,7 +2093,7 @@ imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
2102 2093
 ##	##batch <- unique(object$batch)
2103 2094
 ##	batch <- unique(batch(object))
2104 2095
 ##	tau2A <- getParam(object, "tau2A", batch)
2105
-##	tau2A[index.BB] <- rowMAD(log2(x), log2(x), na.rm=TRUE)^2	
2096
+##	tau2A[index.BB] <- rowMAD(log2(x), log2(x), na.rm=TRUE)^2
2106 2097
 ##	DF <- Ns[, "BB"]-1
2107 2098
 ##	DF[DF < 1] <- 1
2108 2099
 ##	med <- median(tau2A, na.rm=TRUE)
... ...
@@ -2110,17 +2101,17 @@ imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
2110 2101
 ##	tau2A[is.na(tau2A) & isSnp(object)] <- med
2111 2102
 ##	object <- pr(object, "tau2A", batch, tau2A)
2112 2103
 ##
2113
-##	sig2B <- getParam(object, "sig2B", batch)	
2104
+##	sig2B <- getParam(object, "sig2B", batch)
2114 2105
 ##	x <- BB.B[index.BB, ]
2115
-##	sig2B[index.BB] <- rowMAD(log2(x), log2(x), na.rm=TRUE)^2	
2106
+##	sig2B[index.BB] <- rowMAD(log2(x), log2(x), na.rm=TRUE)^2
2116 2107
 ##	med <- median(sig2B, na.rm=TRUE)
2117 2108
 ##	sig2B <- (sig2B * DF  +  med * DF.PRIOR)/(DF.PRIOR + DF)
2118 2109
 ##	sig2B[is.na(sig2B) & isSnp(object)] <- med
2119
-##	object <- pr(object, "sig2B", batch, sig2B)	
2110
+##	object <- pr(object, "sig2B", batch, sig2B)
2120 2111
 ##
2121
-##	tau2B <- getParam(object, "tau2B", batch)	
2112
+##	tau2B <- getParam(object, "tau2B", batch)
2122 2113
 ##	x <- AA.B[index.AA, ]
2123
-##	tau2B[index.AA] <- rowMAD(log2(x), log2(x), na.rm=TRUE)^2		
2114
+##	tau2B[index.AA] <- rowMAD(log2(x), log2(x), na.rm=TRUE)^2
2124 2115
 ##	DF <- Ns[, "AA"]-1
2125 2116
 ##	DF[DF < 1] <- 1
2126 2117
 ##	med <- median(tau2B, na.rm=TRUE)
... ...
@@ -2128,9 +2119,9 @@ imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
2128 2119
 ##	tau2B[is.na(tau2B) & isSnp(object)] <- med
2129 2120
 ##	object <- pr(object, "tau2B", batch, tau2B)
2130 2121
 ##
2131
-##	sig2A <- getParam(object, "sig2A", batch)	
2122
+##	sig2A <- getParam(object, "sig2A", batch)
2132 2123
 ##	x <- AA.A[index.AA, ]
2133
-##	sig2A[index.AA] <- rowMAD(log2(x), log2(x), na.rm=TRUE)^2##var(log(IA)|AA)	
2124
+##	sig2A[index.AA] <- rowMAD(log2(x), log2(x), na.rm=TRUE)^2##var(log(IA)|AA)
2134 2125
 ##	med <- median(sig2A, na.rm=TRUE)
2135 2126
 ##	sig2A <- (sig2A * DF  +  med * DF.PRIOR)/(DF.PRIOR + DF)
2136 2127
 ##	sig2A[is.na(sig2A) & isSnp(object)] <- med
... ...
@@ -2147,7 +2138,7 @@ imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
2147 2138
 ##		med <- median(corr, na.rm=TRUE)
2148 2139
 ##		corr <- (corr*DF  +  med * DF.PRIOR)/(DF.PRIOR + DF)
2149 2140
 ##		corr[is.na(corr) & isSnp(object)] <- med
2150
-##		object <- pr(object, "corr", batch, corr)		
2141
+##		object <- pr(object, "corr", batch, corr)
2151 2142
 ##	}
2152 2143
 ##	corrB.AA <- getParam(object, "corrB.AA", batch)
2153 2144
 ##	backgroundB <- AA.B[index.AA, ]
... ...
@@ -2160,7 +2151,7 @@ imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
2160 2151
 ##	corrB.AA[is.na(corrB.AA) & isSnp(object)] <- med
2161 2152
 ##	object <- pr(object, "corrB.AA", batch, corrB.AA)
2162 2153
 ##
2163
-##	corrA.BB <- getParam(object, "corrA.BB", batch)	
2154
+##	corrA.BB <- getParam(object, "corrA.BB", batch)
2164 2155
 ##	backgroundA <- BB.A[index.BB, ]
2165 2156
 ##	signalB <- BB.B[index.BB, ]
2166 2157
 ##	corrA.BB[index.BB] <- rowCors(backgroundA, signalB, na.rm=TRUE)
... ...
@@ -2195,7 +2186,7 @@ imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
2195 2186
 ##			IB <- muB[, -4]
2196 2187
 ##			vA <- vA[, -4]
2197 2188
 ##			vB <- vB[, -4]
2198
-##			Np <- Ns[, -4] 
2189
+##			Np <- Ns[, -4]
2199 2190
 ##		} else{
2200 2191
 ##			IA <- muA
2201 2192
 ##			IB <- muB
... ...
@@ -2203,7 +2194,7 @@ imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
2203 2194
 ##			vB <- vB
2204 2195
 ##			Np <- Ns
2205 2196
 ##		}
2206
-##		
2197
+##
2207 2198
 ##	}
2208 2199
 ##	Np[Np < 1] <- 1
2209 2200
 ##	vA2 <- vA^2/Np
... ...
@@ -2214,7 +2205,7 @@ imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
2214 2205
 ##	YB <- IB*wB
2215 2206
 ##	##update lm.coefficients stored in object
2216 2207
 ##	object <- nuphiAllele(object, allele="A", Ystar=YA, W=wA, tmp.objects=tmp.objects, cnOptions=cnOptions)
2217
-##	object <- nuphiAllele(object, allele="B", Ystar=YB, W=wB, tmp.objects=tmp.objects, cnOptions=cnOptions)	
2208
+##	object <- nuphiAllele(object, allele="B", Ystar=YB, W=wB, tmp.objects=tmp.objects, cnOptions=cnOptions)
2218 2209
 ##	##---------------------------------------------------------------------------
2219 2210
 ##	##Estimate crosshyb using X chromosome and sequence information
2220 2211
 ##	##---------------------------------------------------------------------------
... ...
@@ -2241,7 +2232,7 @@ imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
2241 2232
 ##	vA <- tmp.objects[["vA"]]
2242 2233
 ##	vB <- tmp.objects[["vB"]]
2243 2234
 ##	Ns <- tmp.objects[["Ns"]]
2244
-##	
2235
+##
2245 2236
 ##	nuA <- getParam(object, "nuA", batch)
2246 2237
 ##	nuB <- getParam(object, "nuB", batch)
2247 2238
 ####	nuA.se <- getParam(object, "nuA.se", batch)
... ...
@@ -2260,12 +2251,12 @@ imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
2260 2251
 ##	##---------------------------------------------------------------------------
2261 2252
 ##	if(CHR == 23){
2262 2253
 ##		phiAX <- getParam(object, "phiAX", batch)  ##nonspecific hybridization coef
2263
-##		phiBX <- getParam(object, "phiBX", batch)  ##nonspecific hybridization coef			
2264
-##		phistar <- phiBX/phiA  
2254
+##		phiBX <- getParam(object, "phiBX", batch)  ##nonspecific hybridization coef
2255
+##		phistar <- phiBX/phiA
2265 2256
 ####		tmp <- (B-nuB - phistar*A + phistar*nuA)/phiB
2266 2257
 ####		copyB <- tmp/(1-phistar*phiAX/phiB)
2267 2258
 ####		copyA <- (A-nuA-phiAX*copyB)/phiA
2268
-####		CB(object) <- copyB  ## multiplies by 100 and converts to integer 
2259
+####		CB(object) <- copyB  ## multiplies by 100 and converts to integer
2269 2260
 ####		CA(object) <- copyA
2270 2261
 ##	} else{
2271 2262
 ####		CA(object) <- matrix((1/phiA*(A-nuA)), nrow(A), ncol(A))
... ...
@@ -2290,21 +2281,21 @@ imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
2290 2281
 ##	sig2A <- getParam(object, "sig2A", batch)[I]
2291 2282
 ##	sig2B <- getParam(object, "sig2B", batch)[I]
2292 2283
 ##	tau2A <- getParam(object, "tau2A", batch)[I]
2293
-##	tau2B <- getParam(object, "tau2B", batch)[I]	
2284
+##	tau2B <- getParam(object, "tau2B", batch)[I]
2294 2285
 ##	corrA.BB <- getParam(object, "corrA.BB", batch)[I]
2295 2286
 ##	corrB.AA <- getParam(object, "corrB.AA", batch)[I]
2296
-##	corr <- getParam(object, "corr", batch)[I]		
2287
+##	corr <- getParam(object, "corr", batch)[I]
2297 2288
 ##	nuA <- getParam(object, "nuA", batch)[I]
2298
-##	nuB <- getParam(object, "nuB", batch)[I]	
2289
+##	nuB <- getParam(object, "nuB", batch)[I]
2299 2290
 ##	phiA <- getParam(object, "phiA", batch)[I]
2300 2291
 ##	phiB <- getParam(object, "phiB", batch)[I]
2301 2292
 ##	normal <- tmp.objects[["normal"]][I, ]
2302 2293
 ##	prior.prob <- cnOptions$prior.prob
2303
-##	emit <- array(NA, dim=c(nrow(A), ncol(A), 10))##SNPs x sample x 'truth'	
2294
+##	emit <- array(NA, dim=c(nrow(A), ncol(A), 10))##SNPs x sample x 'truth'
2304 2295
 ##	lA <- log2(A)
2305
-##	lB <- log2(B)	
2296
+##	lB <- log2(B)
2306 2297
 ##	X <- cbind(lA, lB)
2307
-##	counter <- 1##state counter								
2298
+##	counter <- 1##state counter
2308 2299
 ##	for(CT in 0:3){
2309 2300
 ##		for(CA in 0:CT){
2310 2301
 ##			cat(".")
... ...
@@ -2318,7 +2309,7 @@ imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
2318 2309
 ##			if(CHR == 23){
2319 2310
 ##				##means <- cbind(suppressWarnings(log2(nuA[, p]+CA*phiA[, p] + CB*phiAx[, p])), suppressWarnings(log2(nuB[, p]+CB*phiB[, p] + CA*phiBx[, p])))
2320 2311
 ##				meanA <- suppressWarnings(log2(nuA+CA*phiA + CB*phiAX))
2321
-##				meanB <- suppressWarnings(log2(nuB+CB*phiB + CA*phiBX))				
2312
+##				meanB <- suppressWarnings(log2(nuB+CB*phiB + CA*phiBX))
2322 2313
 ##			} else{
2323 2314
 ##				##means <- cbind(suppressWarnings(log2(nuA+CA*phiA)), suppressWarnings(log2(nuB+CB*phiB)))
2324 2315
 ##				meanA <- suppressWarnings(log2(nuA+CA*phiA))
... ...
@@ -2400,7 +2391,7 @@ imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
2400 2391
 ##		  prior.prob,
2401 2392
 ##		  MIN.SAMPLES,
2402 2393
 ##		  verbose){
2403
-##	
2394
+##
2404 2395
 ##}
2405 2396
 
2406 2397
 ##bias2 <- function(idxBatch,
... ...
@@ -2451,7 +2442,7 @@ imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
2451 2442
 ##
2452 2443
 ##		index <- which(prUp > 0.05 & prUp > prDn)
2453 2444
 ##		##if proportion up greater than 5%, trim the high cn est.
2454
-##		norm[index, k] <- argmax.cn[index, ] > 3 
2445
+##		norm[index, k] <- argmax.cn[index, ] > 3
2455 2446
 ##
2456 2447
 ##		index <- which(prDn > 0.05 & prDn > prUp)
2457 2448
 ##		norm[index, k] <- argmax.cn[index, ] < 3
... ...
@@ -2494,7 +2485,7 @@ imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
2494 2485
 ##		 prior.prob=prior.prob,
2495 2486
 ##		 emit=emit,
2496 2487
 ##		 MIN.SAMPLES=MIN.SAMPLES,
2497
-##		 verbose=verbose)	
2488
+##		 verbose=verbose)
2498 2489
 ##}
2499 2490
 
2500 2491
 
... ...
@@ -2512,7 +2503,7 @@ imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
2512 2503
 ##	nuA <- getParam(object, "nuA", batch)
2513 2504
 ##	phiA <- getParam(object, "phiA", batch)
2514 2505
 ##	prior.prob <- cnOptions$prior.prob
2515
-##	emit <- array(NA, dim=c(nrow(A), ncol(A), 4))##SNPs x sample x 'truth'	
2506
+##	emit <- array(NA, dim=c(nrow(A), ncol(A), 4))##SNPs x sample x 'truth'
2516 2507
 ##	lT <- log2(A)
2517 2508
 ##	I <- isSnp(object)
2518 2509
 ##	counter <- 1 ##state counter
... ...
@@ -2533,7 +2524,7 @@ imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
2533 2524
 ##		counter <- counter+1
2534 2525
 ##	}
2535 2526
 ##	mostLikelyState <- apply(emit, c(1, 2), function(x) order(x, decreasing=TRUE)[1])
2536
-##	
2527
+##
2537 2528
 ##	if(CHR == 23){
2538 2529
 ##		## the state index for male on chromosome 23  is 2
2539 2530
 ##		## add 1 so that the state index is 3 for 'normal' state
... ...
@@ -2563,11 +2554,11 @@ imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
2563 2554
 ##getParams <- function(object, batch){
2564 2555
 ##	##batch <- unique(object$batch)
2565 2556
 ##	batch <- unique(batch(object))
2566
-##	if(length(batch) > 1) stop("batch variable not unique")		
2557
+##	if(length(batch) > 1) stop("batch variable not unique")
2567 2558
 ##	nuA <- as.numeric(fData(object)[, match(paste("nuA", batch, sep="_"), fvarLabels(object))])
2568
-##	nuB <- as.numeric(fData(object)[, match(paste("nuB", batch, sep="_"), fvarLabels(object))])	
2559
+##	nuB <- as.numeric(fData(object)[, match(paste("nuB", batch, sep="_"), fvarLabels(object))])
2569 2560
 ##	phiA <- as.numeric(fData(object)[, match(paste("phiA", batch, sep="_"), fvarLabels(object))])
2570
-##	phiB <- as.numeric(fData(object)[, match(paste("phiB", batch, sep="_"), fvarLabels(object))])	
2561
+##	phiB <- as.numeric(fData(object)[, match(paste("phiB", batch, sep="_"), fvarLabels(object))])
2571 2562
 ##	tau2A <- as.numeric(fData(object)[, match(paste("tau2A", batch, sep="_"), fvarLabels(object))])
2572 2563
 ##	tau2B <- as.numeric(fData(object)[, match(paste("tau2B", batch, sep="_"), fvarLabels(object))])
2573 2564
 ##	sig2A <- as.numeric(fData(object)[, match(paste("sig2A", batch, sep="_"), fvarLabels(object))])
... ...
@@ -2609,17 +2600,17 @@ imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
2609 2600
 ##	phiB <- getParam(object, "phiB", batch)
2610 2601
 ##	if(!all(is.na(phiB))){
2611 2602
 ##		phiB[phiB < MIN.PHI] <- MIN.PHI
2612
-##		object <- pr(object, "phiB", batch, phiB)		
2603
+##		object <- pr(object, "phiB", batch, phiB)
2613 2604
 ##	}
2614
-##	phiAX <- as.numeric(getParam(object, "phiAX", batch))	
2605
+##	phiAX <- as.numeric(getParam(object, "phiAX", batch))
2615 2606
 ##	if(!all(is.na(phiAX))){
2616 2607
 ##		phiAX[phiAX < MIN.PHI] <- MIN.PHI
2617
-##		object <- pr(object, "phiAX", batch, phiAX)		
2608
+##		object <- pr(object, "phiAX", batch, phiAX)
2618 2609
 ##	}
2619 2610
 ##	phiBX <- as.numeric(getParam(object, "phiBX", batch))
2620 2611
 ##	if(!all(is.na(phiBX))){
2621 2612
 ##		phiBX[phiBX < MIN.PHI] <- MIN.PHI
2622
-##		object <- pr(object, "phiBX", batch, phiBX)			
2613
+##		object <- pr(object, "phiBX", batch, phiBX)
2623 2614
 ##	}
2624 2615
 ##	return(object)
2625 2616
 ##}
... ...
@@ -2655,11 +2646,11 @@ imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
2655 2646
 ##		verbose <- cnOptions$verbose
2656 2647
 ##		##if(verbose) message("Thresholding nu and phi")
2657 2648
 ##		object <- thresholdModelParams(object, cnOptions)
2658
-##	}		
2659
-##	##if(verbose) message("\nAllele specific copy number")	
2649
+##	}
2650
+##	##if(verbose) message("\nAllele specific copy number")
2660 2651
 ##	object <- polymorphic(object, cnOptions, tmp.objects)
2661 2652
 ##	if(any(!isSnp(object))){  ##there are nonpolymorphic probes
2662
-##		##if(verbose) message("\nCopy number for nonpolymorphic probes...")	
2653
+##		##if(verbose) message("\nCopy number for nonpolymorphic probes...")
2663 2654
 ##		object <- nonpolymorphic(object, cnOptions, tmp.objects)
2664 2655
 ##	}
2665 2656
 ##	##---------------------------------------------------------------------------
... ...
@@ -2677,9 +2668,9 @@ imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
2677 2668
 ##	object <- pr(object, "phiB", PLATE, getParam(object, "phiB", PLATE))
2678 2669
 ##	object <- pr(object, "phiB.se", PLATE, getParam(object, "phiB.se", PLATE))
2679 2670
 ##	object <- pr(object, "tau2A", PLATE, getParam(object, "tau2A", PLATE))
2680
-##	object <- pr(object, "tau2B", PLATE, getParam(object, "tau2B", PLATE))				
2671
+##	object <- pr(object, "tau2B", PLATE, getParam(object, "tau2B", PLATE))
2681 2672
 ##	object <- pr(object, "sig2A", PLATE, getParam(object, "sig2A", PLATE))
2682
-##	object <- pr(object, "sig2B", PLATE, getParam(object, "sig2B", PLATE))		
2673
+##	object <- pr(object, "sig2B", PLATE, getParam(object, "sig2B", PLATE))
2683 2674
 ##	object <- pr(object, "phiAX", PLATE, as.numeric(getParam(object, "phiAX", PLATE)))
2684 2675
 ##	object <- pr(object, "phiBX", PLATE, as.numeric(getParam(object, "phiBX", PLATE)))
2685 2676
 ##	object <- pr(object, "corr", PLATE, getParam(object, "corr", PLATE))
... ...
@@ -2696,7 +2687,7 @@ imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
2696 2687
 
2697 2688
 ## constructors for Illumina platform
2698 2689
 constructIlluminaFeatureData <- function(gns, cdfName){
2699
-	pkgname <- paste(cdfName, "Crlmm", sep="")	
2690
+	pkgname <- paste(cdfName, "Crlmm", sep="")
2700 2691
 	path <- system.file("extdata", package=pkgname)
2701 2692
 	load(file.path(path, "cnProbes.rda"))
2702 2693
 	load(file.path(path, "snpProbes.rda"))
... ...
@@ -2712,7 +2703,7 @@ constructIlluminaFeatureData <- function(gns, cdfName){
2712 2703
 	colnames(mapping) <- c("chromosome", "position", "isSnp")
2713 2704
 	new("AnnotatedDataFrame",
2714 2705
 	    data=data.frame(mapping),
2715
-	    varMetadata=data.frame(labelDescription=colnames(mapping)))	
2706
+	    varMetadata=data.frame(labelDescription=colnames(mapping)))
2716 2707
 }
2717 2708
 constructIlluminaAssayData <- function(np, snp, object, storage.mode="environment", order.index){
2718 2709
 	stopifnot(identical(snp$gns, featureNames(object)))
... ...
@@ -2746,14 +2737,14 @@ constructIlluminaCNSet <- function(crlmmResult,
2746 2737
 	load(file.path(path, "snpFile.rda"))
2747 2738
 	res <- get("res")
2748 2739
 	load(file.path(path, "cnFile.rda"))
2749
-	cnAB <- get("cnAB")	
2740
+	cnAB <- get("cnAB")
2750 2741
 	fD <- constructIlluminaFeatureData(c(res$gns, cnAB$gns), cdfName="human370v1c")
2751 2742
 	new.order <- order(fD$chromosome, fD$position)
2752 2743
 	fD <- fD[new.order, ]
2753 2744
 	aD <- constructIlluminaAssayData(cnAB, res, crlmmResult, order.index=new.order)
2754 2745
 	##protocolData(crlmmResult)$batch <- vector("integer", ncol(crlmmResult))
2755 2746
 	batch <- vector("integer", ncol(crlmmResult))
2756
-	container <- new("CNSet", 
2747
+	container <- new("CNSet",
2757 2748
 			 call=aD[["call"]],
2758 2749
 			 callProbability=aD[["callProbability"]],
2759 2750
 			 alleleA=aD[["alleleA"]],
... ...
@@ -2764,11 +2755,11 @@ constructIlluminaCNSet <- function(crlmmResult,
2764 2755
 			 batch=batch,
2765 2756
 			 annotation="human370v1c")
2766 2757
 ##	lM(container) <- initializeParamObject(list(featureNames(container), unique(protocolData(container)$batch)))
2767
-	lM(container) <- initializeLmFrom(container)
2758
+##	lM(container) <- initializeLmFrom(container)
2768 2759
 	container
2769 2760
 }
2770
-				   
2771
-				   
2761
+
2762
+
2772 2763
 
2773 2764
 ellipseCenters <- function(object, index, allele, batch, log.it=TRUE){
2774 2765
 	ubatch <- unique(batch(object))[batch]
... ...
@@ -2788,7 +2779,7 @@ ellipseCenters <- function(object, index, allele, batch, log.it=TRUE){
2788 2779
 	names(centers) <- nms
2789 2780
 	fns <- featureNames(object)[index]
2790 2781
 	centers$fns <- fns
2791
-	return(centers)	
2782
+	return(centers)
2792 2783
 }
2793 2784
 
2794 2785
 computeCN <- function(object,
... ...
@@ -2831,12 +2822,12 @@ computeCN <- function(object,
2831 2822
 			save(flags, file=file.path(ldPath(), "flags.snps.rda"))
2832 2823
 		} else{
2833 2824
 			load(file.path(ldPath(), "flags.snps.rda"))
2834
-		} 
2825
+		}
2835 2826
 		if(verbose) message("Estimating allele-specific copy number at autosomal SNPs")
2836 2827
 		FUN <- fit.lm1
2837 2828
 	}
2838 2829
 	if(type=="autosome.nps"){
2839
-		marker.index <- which(chromosome(object) < 23 & !isSnp(object) & !is.na(chromosome(object)))		
2830
+		marker.index <- which(chromosome(object) < 23 & !isSnp(object) & !is.na(chromosome(object)))
2840 2831
 		##marker.index <- which(chromosome(object) < 23 & isSnp(object) & !is.na(chromosome(object)))
2841 2832
 		nr <- length(marker.index)
2842 2833
 		Ns <- initializeBigMatrix("Ns", nr, 5)
... ...
@@ -2852,7 +2843,7 @@ computeCN <- function(object,
2852 2843
 			save(flags, file=file.path(ldPath(), "flags.nps.rda"))
2853 2844
 		} else{
2854 2845
 			load(file.path(ldPath(), "flags.nps.rda"))
2855
-		} 
2846
+		}
2856 2847
 		if(verbose) message("Estimating allele-specific copy number at autosomal SNPs")
2857 2848
 		FUN <- fit.lm2
2858 2849
 	}
... ...
@@ -2873,12 +2864,12 @@ computeCN <- function(object,
2873 2864
 			save(flags, file=file.path(ldPath(), "flags.X.snps.rda"))
2874 2865
 		} else{
2875 2866
 			load(file.path(ldPath(), "flags.X.snps.rda"))
2876
-		} 
2867
+		}
2877 2868
 		if(verbose) message("Estimating allele-specific copy number at autosomal SNPs")
2878 2869
 		FUN <- fit.lm3
2879 2870
 	}
2880 2871
 	if(type=="X.nps"){
2881
-		marker.index <- which(chromosome(object) == 23 & !isSnp(object) & !is.na(chromosome(object)))		
2872
+		marker.index <- which(chromosome(object) == 23 & !isSnp(object) & !is.na(chromosome(object)))
2882 2873
 		nr <- length(marker.index)
2883 2874
 		Ns <- initializeBigMatrix("Ns", nr, 5)
2884 2875
 		colnames(Ns) <- c("A", "B", "AA", "AB", "BB")
... ...
@@ -2893,7 +2884,7 @@ computeCN <- function(object,
2893 2884
 			save(flags, file=file.path(ldPath(), "flags.X.nps.rda"))
2894 2885
 		} else{
2895 2886
 			load(file.path(ldPath(), "flags.X.nps.rda"))
2896
-		} 
2887
+		}
2897 2888
 		if(verbose) message("Estimating allele-specific copy number at autosomal SNPs")
2898 2889
 		FUN <- fit.lm4
2899 2890
 	}
... ...
@@ -17,8 +17,6 @@ linearParamElementReplace <- function(obj, elt, value) {
17 17
 }
18 18
 
19 19
 
20
-
21
-
22 20
 setMethod("nuA", signature=signature(object="CNSet"), function(object) nu(object, "A"))
23 21
 setMethod("nuB", signature=signature(object="CNSet"), function(object) nu(object, "B"))
24 22
 setMethod("phiA", signature=signature(object="CNSet"), function(object) phi(object, "A"))
... ...
@@ -31,122 +29,79 @@ setMethod("corrAA", signature=signature(object="CNSet"), function(object) corr(o
31 29
 setMethod("corrAB", signature=signature(object="CNSet"), function(object) corr(object, "AB"))
32 30
 setMethod("corrBB", signature=signature(object="CNSet"), function(object) corr(object, "BB"))
33 31
 
34
-setReplaceMethod("nuA", signature=signature(object="CNSet", value="matrix"), 
32
+setReplaceMethod("nuA", signature=signature(object="CNSet", value="matrix"),
35 33
 	  function(object, value){
36 34
 		  linearParamElementReplace(object, "nuA", value)
37 35
 	  })
38 36
 
39
-setReplaceMethod("nuB", signature=signature(object="CNSet", value="matrix"), 
37
+setReplaceMethod("nuB", signature=signature(object="CNSet", value="matrix"),
40 38
 	  function(object, value){
41
-		  linearParamElementReplace(object, "nuB", value)		  
39
+		  linearParamElementReplace(object, "nuB", value)
42 40
 })
43 41
 
44
-setReplaceMethod("phiA", signature=signature(object="CNSet", value="matrix"), 
42
+setReplaceMethod("phiA", signature=signature(object="CNSet", value="matrix"),
45 43
 	  function(object, value){
46
-		  linearParamElementReplace(object, "phiA", value)		  
44
+		  linearParamElementReplace(object, "phiA", value)
47 45
 })
48 46
 
49
-setReplaceMethod("phiB", signature=signature(object="CNSet", value="matrix"), 
47
+setReplaceMethod("phiB", signature=signature(object="CNSet", value="matrix"),
50 48
 	  function(object, value){
51
-		  linearParamElementReplace(object, "phiB", value)		  
49
+		  linearParamElementReplace(object, "phiB", value)
52 50
 })
53 51
 
54
-setReplaceMethod("sigma2A", signature=signature(object="CNSet", value="matrix"), 
52
+setReplaceMethod("sigma2A", signature=signature(object="CNSet", value="matrix"),
55 53
 	  function(object, value){
56
-		  linearParamElementReplace(object, "sig2A", value)		  
54
+		  linearParamElementReplace(object, "sig2A", value)
57 55
 })
58 56
 
59
-setReplaceMethod("sigma2B", signature=signature(object="CNSet", value="matrix"), 
57
+setReplaceMethod("sigma2B", signature=signature(object="CNSet", value="matrix"),
60 58
 	  function(object, value){
61
-		  linearParamElementReplace(object, "sig2B", value)		  
59
+		  linearParamElementReplace(object, "sig2B", value)
62 60
 })
63 61
 
64
-setReplaceMethod("tau2A", signature=signature(object="CNSet", value="matrix"), 
62
+setReplaceMethod("tau2A", signature=signature(object="CNSet", value="matrix"),
65 63
 	  function(object, value){
66
-		  linearParamElementReplace(object, "tau2A", value)		  
64
+		  linearParamElementReplace(object, "tau2A", value)
67 65
 })
68 66
 
69
-setReplaceMethod("tau2B", signature=signature(object="CNSet", value="matrix"), 
67
+setReplaceMethod("tau2B", signature=signature(object="CNSet", value="matrix"),
70 68
 	  function(object, value){
71
-		  linearParamElementReplace(object, "tau2B", value)		  
69
+		  linearParamElementReplace(object, "tau2B", value)
72 70
 })
73 71
 
74
-setReplaceMethod("corrAA", signature=signature(object="CNSet", value="matrix"), 
72
+setReplaceMethod("corrAA", signature=signature(object="CNSet", value="matrix"),
75 73
 	  function(object, value){
76
-		  linearParamElementReplace(object, "corrAA", value)		  
74
+		  linearParamElementReplace(object, "corrAA", value)
77 75
 })
78 76
 
79
-setReplaceMethod("corrAB", signature=signature(object="CNSet", value="matrix"), 
77
+setReplaceMethod("corrAB", signature=signature(object="CNSet", value="matrix"),
80 78
 	  function(object, value){
81
-		  linearParamElementReplace(object, "corrAB", value)		  
79
+		  linearParamElementReplace(object, "corrAB", value)
82 80
 })
83 81
 
84
-setReplaceMethod("corrBB", signature=signature(object="CNSet", value="matrix"), 
82
+setReplaceMethod("corrBB", signature=signature(object="CNSet", value="matrix"),
85 83
 	  function(object, value){
86
-		  linearParamElementReplace(object, "corrBB", value)		  
84
+		  linearParamElementReplace(object, "corrBB", value)
87 85
 })
88 86
 
89 87
 
90
-##setValidity("CNSet",
91
-##	    function(object){
92
-##		    if(!"batch" %in% varLabels(protocolData(object)))
93
-##			    return("'batch' not defined in protocolData")
94
-##		    if(!"chromosome" %in% fvarLabels(object))
95
-##			    return("'chromosome' not defined in featureData")
96
-##		    if(!"position" %in% fvarLabels(object))
97
-##			    return("'position' not defined in featureData")
98
-##		    if(!"isSnp" %in% fvarLabels(object))
99
-##			    return("'isSnp' not defined in featureData")
100
-##		    return(TRUE)
101
-##	    })
102
-
103
-setMethod("totalCopyNumber", "CNSet", function(object, i, j){
104
-	if(missing(i) & missing(j)){
105
-		if(inherits(CA(object), "ff") | inherits(CA(object), "ffdf")) stop("Must specify i and/or j for ff objects")
106
-	}
107
-	if(missing(i) & !missing(j)){
108
-		snp.index <- which(isSnp(object))	
109
-		cn.total <- as.matrix(CA(cnSet)[, j])
110
-		cb <- as.matrix(CB(cnSet)[snp.index, j]	)
111
-		cn.total[snp.index, ] <- cn.total[snp.index, ] + cb		
112
-	}
113
-	if(!missing(i) & missing(j)){
114
-		snp.index <- intersect(which(isSnp(object)), i)
115
-		cn.total <- as.matrix(CA(cnSet)[i, ])
116
-		cb <- as.matrix(CB(cnSet)[snp.index, ])	
117
-		cn.total[snp.index, ] <- cn.total[snp.index, ] + cb				
118
-	}
119
-	if(!missing(i) & !missing(j)){
120
-		snp.index <- intersect(which(isSnp(object)), i)		
121
-		cn.total <- as.matrix(CA(cnSet)[i, j])	
122
-		cb <- as.matrix(CB(cnSet)[snp.index, j])
123
-		cn.total[snp.index, ] <- cn.total[snp.index, ] + cb
124
-	}
125
-	cn.total <- cn.total/100
126
-	dimnames(cn.total) <- NULL
127
-	return(cn.total)
128
-})
129
-
130 88
 ##setMethod("ellipse", "CNSet", function(x, copynumber, batch, ...){
131 89
 ##	ellipse.CNSet(x, copynumber, batch, ...)
132 90
 ##})
133 91
 
134 92
 
135
-
136 93
 ACN <- function(object, allele, i , j){
137
-	if(missing(i) & missing(j)){
138
-		if(inherits(A(object), "ff") | inherits(A(object), "ffdf")) stop("Must specify i and/or j for ff objects")
139
-	}
94
+	bns <- batchNames(object)
95
+	acn <- list()
140 96
 	if(missing(i) & !missing(j)){
141 97
 		## calculate ca only for batches indexed by j
142
-		ubatch <- unique(batch(object))
143 98
 		batches <- unique(batch(object)[j])
144 99
 		for(k in seq_along(batches)){
145
-			l <- match(batches[k], ubatch)
100
+			l <- match(batches[k], bns)
146 101
 			bg <- nu(object, allele)[, l]
147
-			sl <- phi(object, allele)[, l]
102
+			slope <- phi(object, allele)[, l]
148 103
 			I <- allele(object, allele)[, j]
149
-			acn <- 1/sl*(I - bg)				  
104
+			acn[[k]] <- 1/slope*(I - bg)
150 105
 		}
151 106
 	}
152 107
 	if(!missing(i) & missing(j)){
... ...
@@ -155,75 +110,93 @@ ACN <- function(object, allele, i , j){
155 110
 		for(k in seq_along(batches)){
156 111
 			##bb <- batches[k]
157 112
 			bg <- nu(object, allele)[i, k]
158
-			sl <- phi(object, allele)[i, k]
113
+			slope <- phi(object, allele)[i, k]
159 114
 			I <- allele(object, allele)[i, j]
160
-			acn <- 1/sl*(I - bg)
115
+			acn[[k]] <- 1/slope*(I - bg)
161 116
 		}
162 117
 	}
163 118
 	if(!missing(i) & !missing(j)){
164 119
 		ubatch <- unique(batch(object))
165 120
 		batches <- unique(batch(object)[j])
121
+		acn <- list()
166 122
 		for(k in seq_along(batches)){
167 123
 			l <- match(batches[k], ubatch)
168 124
 			bg <- nu(object, allele)[i, l]
169
-			sl <- phi(object, allele)[i, l]
125
+			slope <- phi(object, allele)[i, l]
170 126
 			I <- allele(object, allele)[i, j]
171
-			acn <- 1/bg*(I - sl)				  
172
-		}			  
127
+			acn[[k]] <- 1/slope*(I - bg)
128
+		}
173 129
 	}
130
+	if(length(acn) > 1) acn <- do.call("cbind", acn)
131
+	if(length(acn) == 1) acn <- acn[[1]]
174 132
 	return(acn)
175 133
 }
176 134
 
177
-setMethod("CA",
178
-	  signature=signature(object="CNSet", i="integerOrMissing", j="integerOrMissing"),
179
-	  function(object, i, j) {
180
-		  ca <- ACN(object, allele="A", i, j)
135
+setMethod("CA", signature=signature(object="CNSet"),
136
+	  function(object, ...){
137
+		  ca <- ACN(object, allele="A", ...)
181 138
 		  return(ca)
182 139
 	  })
183
-setMethod("CB",
184
-	  signature=signature(object="CNSet", i="integerOrMissing", j="integerOrMissing"),
185
-	  function(object, i, j) {
186
-		  cb <- ACN(object, allele="B", i, j)
140
+setMethod("CB", signature=signature(object="CNSet"),
141
+	  function(object, ...) {
142
+		  cb <- ACN(object, allele="B", ...)
187 143
 		  return(cb)
188 144
 	  })
189 145
 
190
-setMethod("totalCopyNumber",
191
-	  signature=signature(object="CNSet", i="integerOrMissing", j="integerOrMissing"),
192
-	  function(object, i, j, ...){
193
-	if(missing(i) & missing(j)){
194
-		if(inherits(A(object), "ff") | inherits(A(object), "ffdf")) stop("Must specify i and/or j for ff objects")
195
-	}
196
-	if(missing(i) & !missing(j)){
197
-		snp.index <- which(isSnp(object))	
198
-		cn.total <- as.matrix(CA(object)[, j])
199
-		if(length(snp.index) > 0){
200
-			cb <- as.matrix(CB(object)[snp.index, j])
201
-			snps <- (1:nrow(cn.total))[i %in% snp.index]
202
-			cn.total[snps, ] <- cn.total[snps, j] + cb				
203
-		}
204
-	}
205
-	if(!missing(i) & missing(j)){
206
-		snp.index <- intersect(which(isSnp(object)), i)
207
-		cn.total <- as.matrix(CA(object)[i, ])
208
-		if(length(snp.index) > 0){
209
-			cb <- as.matrix(CB(object)[snp.index, ])
210
-			snps <- (1:nrow(cn.total))[i %in% snp.index]
211
-			cn.total[snps, ] <- cn.total[snps, ] + cb				
212
-		}
146
+##setMethod("totalCopyNumber", signature=signature(object="CNSet"),
147
+totalCopyNumber <- function(object, ..., verbose=TRUE, dimnames=FALSE){
148
+	ca <- CA(object, ...)
149
+	cb <- CB(object, ...)
150
+	is.snp <- isSnp(object)
151
+	dotArgs <- list(...)
152
+	if("i" %in% names(dotArgs)){
153
+		i <- dotArgs[["i"]]
154
+		np.index <- which(!is.snp[i])
155
+		if(length(np.index) > 0) cb[np.index, ] <- 0
156
+	} else {
157
+		np.index <- which(!is.snp)
158
+		if(length(np.index) > 0) cb[np.index, ] <- 0
213 159
 	}
214
-	if(!missing(i) & !missing(j)){
215
-		snp.index <- intersect(which(isSnp(object)), i)		
216
-		cn.total <- as.matrix(CA(object)[i, j])
217
-		if(length(snp.index) > 0){
218
-			cb <- as.matrix(CB(object)[snp.index, j])
219
-			snps <- (1:nrow(cn.total))[i %in% snp.index]
220
-			cn.total[snps, ] <- cn.total[snps, ] + cb
221
-		}
222
-	}
223
-	##cn.total <- cn.total/100
224
-	dimnames(cn.total) <- NULL
225
-	return(cn.total)
226
-})
160
+	return(ca+cb)
161
+}
162
+
163
+##
164
+##
165
+##
166
+##
167
+##		  if(missing(i) & !missing(j)){
168
+##			  ca <- CA(object, i, j)
169
+##			  snp.index <- which(isSnp(object))
170
+##			  cn.total <- as.matrix(CA(object)[, j])
171
+##			  if(length(snp.index) > 0){
172
+##				  cb <- as.matrix(CB(object)[snp.index, j])
173
+##				  snps <- (1:nrow(cn.total))[i %in% snp.index]
174
+##				  cn.total[snps, ] <- cn.total[snps, j] + cb
175
+##			  }
176
+##		  }
177
+##		  if(!missing(i) & missing(j)){
178
+##			  snp.index <- intersect(which(isSnp(object)), i)
179
+##			  cn.total <- as.matrix(CA(object)[i, ])
180
+##			  if(length(snp.index) > 0){
181
+##				  cb <- as.matrix(CB(object)[snp.index, ])
182
+##				  snps <- (1:nrow(cn.total))[i %in% snp.index]
183
+##				  cn.total[snps, ] <- cn.total[snps, ] + cb
184
+##			  }
185
+##		  }
186
+##		  if(!missing(i) & !missing(j)){
187
+##			  snp.index <- intersect(which(isSnp(object)), i)
188
+##			  cn.total <- as.matrix(CA(object)[i, j])
189
+##			  if(length(snp.index) > 0){
190
+##				  cb <- as.matrix(CB(object)[snp.index, j])
191
+##				  snps <- (1:nrow(cn.total))[i %in% snp.index]
192
+##				  cn.total[snps, ] <- cn.total[snps, ] + cb
193
+##			  }
194
+##		  }
195
+##		  ##cn.total <- cn.total/100
196
+##		  dimnames(cn.total) <- NULL
197
+##		  return(cn.total)
198
+##	  })
199
+
227 200
 
228 201
 setReplaceMethod("snpCall", c("CNSet", "ff_or_matrix"),
229 202
                  function(object, ..., value){
... ...
@@ -1,19 +1,19 @@
1
-setMethod("lines", c("CNSetLM"), function(x, y, batch, copynumber, ...){
2
-	linesCNSet(x, y, batch, copynumber, ...)
3
-})
1
+setMethod("lines", signature=signature(x="CNSet"),
2
+	  function(x, y, batch, copynumber, ...){
3
+		  linesCNSet(x, y, batch, copynumber, ...)
4
+	  })
5
+
4 6
 linesCNSet <- function(x, y, batch, copynumber, x.axis="A", ...){
5 7
 	require(ellipse)
6 8
 	object <- x
7 9
 	I <- y
8
-	stopifnot(length(I) == 1) 
9
-	calls.class <- class(calls(object)[[1]])
10
-	ffIsLoaded <- calls.class[1] == "ff_matrix" | calls.class[1] == "ffdf" | calls.class[1]=="ff"
11
-	column <- grep(batch, unique(batch(object)))
10
+	stopifnot(length(I) == 1)
11
+	column <- grep(batch, batchNames(object))
12 12
 	stopifnot(length(column) == 1)
13 13
 	nuA <- nu(object, "A")[I, column]
14 14
 	nuB <- nu(object, "B")[I, column]
15 15
 	phiA <- phi(object, "A")[I, column]
16
-	phiB <- phi(object, "B")[I, column]		
16
+	phiB <- phi(object, "B")[I, column]
17 17
 	tau2A <- tau2(object, "A")[I, column]
18 18
 	tau2B <- tau2(object, "B")[I, column]
19 19
 	sigma2A <- sigma2(object, "A")[I, column]
... ...
@@ -45,5 +45,5 @@ linesCNSet <- function(x, y, batch, copynumber, x.axis="A", ...){
45 45
 					      scale=rev(scale)), ...)
46 46
 			}
47 47
 		}
48
-	}	
48
+	}
49 49
 }
50 50
new file mode 100644
... ...
@@ -0,0 +1,36 @@
1
+\name{CNSet-methods}
2
+\Rdversion{1.1}
3
+\docType{methods}
4
+\alias{CA,CNSet-method}
5
+\alias{CB,CNSet-method}
6
+\alias{lines,CNSet-method}
7
+\alias{totalCopyNumber,CNSet-method}
8
+
9
+\title{crlmm methods for class "CNSet"}
10
+\description{
11
+
12
+	CNSet is a container defined in the oligoClasses package for
13
+	 storing normalized intensities for genotyping platforms,
14
+	 genotype calls, and parameters estimated for copy
15
+	 number. Accessors for data that an object of this class
16
+	 contains are largely defined in the package
17
+	 oligoClasses. CNSet methods that involve more complex
18
+	 calculations that are specific to the crlmm package, such as
19
+	 computing allele-specific copy number, are included in crlmm
20
+	 and described here.
21
+
22
+}
23
+
24
+\section{Methods}{
25
+  \describe{
26
+    \item{CA}{\code{signature(object="CNSet")}: ...}
27
+    \item{CB}{\code{signature(object="CNSet")}: ...}
28
+    \item{lines}{\code{signature(x="CNSet")}: ...}
29
+    \item{totalCopyNumber}{\code{signature(object="CNSet")}: ...}
30
+	 }
31
+}
32
+
33
+\seealso{
34
+  \code{\link{CNSet-class}}, \code{\link{copynumberAccessors}}
35
+}
36
+\keyword{methods}
0 37
deleted file mode 100644
... ...
@@ -1,45 +0,0 @@
1
-\name{batch}
2
-\alias{batch}
3
-\alias{batch,eSet-method}
4
-\title{
5
-	Function to extract batch information.
6
-}
7
-\description{
8
-	Checks the phenoData and protocolData for a variable named
9
-	batch and, if present, returns the vector.
10
-}
11
-\usage{
12
-batch(object)
13
-}
14
-\arguments{
15
-  \item{object}{
16
-  An object extending the \code{eSet} class.
17
-}
18
-}
19
-\details{
20
-
21
-	For copy number estimation, a batch variable must be
22
-	specified.  Currently, we suggest storing this variable in the
23
-	protocolData.
24
-
25
-	Batch represents groups of samples that were processed (DNA
26
-	preparation and collection, PCR amplification, scan date) at
27
-	similar times. Often, the 96 well chemistry plate or scan date
28
-	is a useful surrogate for batch.
29
-
30
-}
31
-\value{
32
-
33
-	Vector indicating batch.
34
-
35
-}
36
-
37
-\author{
38
-R. Scharpf
39
-}
40
-
41
-\seealso{
42
-	\code{\link{genotype}}, \code{\link{genotype2}}
43
-}
44
-\keyword{manip}
45
-
46 0
new file mode 100644
... ...
@@ -0,0 +1,99 @@
1
+\name{copynumberAccessors}
2
+\alias{CA}
3
+\alias{CB}
4
+\alias{totalCopynumber}
5
+
6
+\title{
7
+
8
+	Accessors for allele-specific or total copy number
9
+
10
+}
11
+
12
+
13
+\description{
14
+
15
+These methods can be applied after an object of class \code{CNSet} has
16
+been generated by the \code{crlmmCopynumber} function.
17
+
18
+}
19
+\usage{
20
+CA(object, ...)
21
+CB(object, ...)
22
+totalCopynumber(object,...)
23
+}
24
+
25
+\arguments{
26
+  \item{object}{
27
+  An object of class \code{CNSet}.
28
+}
29
+  \item{\dots}{
30
+
31
+  An additional argument named 'i' can be passed to subset the markers
32
+  and an argument 'j' can be passed to subset the samples.  Other
33
+  arguments are ignored.
34
+
35
+}
36
+}
37
+
38
+\details{
39
+
40
+	These functions can be used to tranlate the normalized
41
+	intensities to the copy number scale.  Plotting the copy
42
+	number estimates as a function of physical position can be
43
+	used to guide downstream algorithms that smooth, as well as to
44
+	assess possible mosaicism.
45
+
46
+}
47
+\value{
48
+	A matrix of allele-specific or total copy number estimates.
49
+}
50
+
51
+\seealso{
52
+	\code{\link{crlmmCopynumber}}, \code{\link{CNSet-class}}, \code{\link{CNSet-methods}}
53
+}
54
+\examples{
55
+data(sample.CNSetLM)
56
+## update to class CNSet
57
+cnSet <- as(sample.CNSetLM, "CNSet")
58
+all(isCurrent(cnSet)) ## is the cnSet object current?
59
+
60
+## --------------------------------------------------
61
+## calculating allele-specific copy number
62
+## --------------------------------------------------
63
+## copy number for allele A, first 5 markers, first 2 samples
64
+(ca <- CA(cnSet, i=1:5, j=1:2))
65
+## copy number for allele B, first 5 markers, first 2 samples
66
+(cb <- CB(cnSet, i=1:5, j=1:2))
67
+## total copy number for first 5 markers, first 2 samples
68
+(cn1 <- ca+cb)
69
+
70
+## total copy number at first 5 nonpolymorphic loci
71
+index <- which(!isSnp(cnSet))[1:5]
72
+cn2 <- CA(cnSet, i=index, j=1:2)
73
+## note, cb is NA at nonpolymorphic loci
74
+(cb <- CB(cnSet, i=index, j=1:2))
75
+## note, ca+cb will give NAs at nonpolymorphic loci
76
+CA(cnSet, i=index, j=1:2) + cb
77
+## A shortcut for total copy number
78
+cn3 <- totalCopyNumber(cnSet, i=1:5, j=1:2)
79
+all.equal(cn3, cn1)
80
+cn4 <- totalCopyNumber(cnSet, i=index, j=1:2)
81
+all.equal(cn4, cn2)
82
+
83
+## markers 1-5, all samples
84
+cn5 <- totalCopyNumber(cnSet, i=1:5)
85
+## all markers, samples 1-5
86
+cn6 <- totalCopyNumber(cnSet, j=1:5)
87
+
88
+## NOTE: subsetting the object before extracting copy number
89
+##       can be very inefficient when the data set is very large,
90
+##       particularly if using ff objects.  IN particular, subsetting
91
+##       the CNSet object will subset all of the assay data elements
92
+##       and all of the elements in the LinearModelParameter slot
93
+\dontrun{
94
+        ## do not do the following
95
+	cn <- CA(cnSet[1:5, ], "A")
96
+}
97
+}
98
+\keyword{manip}
99
+
... ...
@@ -32,12 +32,7 @@ thresholdCopynumber = TRUE)
32 32
 
33 33
 }
34 34
 
35
- \item{chromosome}{Numeric vector indicating which chromosomes to
36
- process (length <= 23). For chromosome X, use 23. A copy number
37
- method for chromosome Y is not yet available.
38
-}
39
-
40
-  \item{MIN.SAMPLES}{ 'Integer'.  The minimum number of samples in a
35
+ \item{MIN.SAMPLES}{ 'Integer'.  The minimum number of samples in a
41 36
   batch.  Bathes with fewer than MIN.SAMPLES are skipped.  Therefore,
42 37
   samples in batches with fewer than MIN.SAMPLES have NA's for the
43 38
   allele-specific copy number and NA's for the linear model
... ...
@@ -2,7 +2,7 @@
2 2
 \alias{sample.CNSetLM}
3 3
 \docType{data}
4 4
 \title{
5
-	Dataset of class 'CNSetLM'	
5
+	Dataset of class 'CNSetLM'
6 6
 }
7 7
 \description{
8 8
 
... ...
@@ -15,20 +15,115 @@
15 15
 
16 16
 	This class has been deprecated.  See example below for how to
17 17
 	update an existing 'CNSetLM' object to class 'CNSet'.
18
- 
18
+
19 19
 	The data illustrates the \code{CNSetLM-class}, with
20 20
 	\code{assayData} containing the quantile-normalized
21 21
 	intensities for the A and B alleles, genotype calls and
22 22
 	confidence scores (call and callProbability), and
23 23
 	allele-specific copy number (CA and CB). The parameters from
24 24
 	the linear model are stored in the lM slot.
25
-	
25
+
26 26
 }
27 27
 \examples{
28 28
 ## class CNSetLM has been deprecated
29 29
 data(sample.CNSetLM)
30 30
 ## update to class CNSet
31 31
 cnSet <- as(sample.CNSetLM, "CNSet")
32
-all(isCurrent(cnSet))
32
+all(isCurrent(cnSet)) ## is the cnSet object current?
33
+
34
+## information on the features
35
+fvarLabels(cnSet)
36
+##
37
+## --------------------------------------------------
38
+## accessors for the feature-level info
39
+## --------------------------------------------------
40
+chromosome(cnSet)[1:5]
41
+position(cnSet)[1:5]
42
+isSnp(cnSet)[1:5]
43
+## 980 nonpolymorphic markers and 1139 polymoprhic markers
44
+table(isSnp(cnSet))
45
+## --------------------------------------------------
46
+
47
+## --------------------------------------------------
48
+## sample-level statistics computed by crlmm
49
+## --------------------------------------------------
50
+varLabels(cnSet)
51
+## accessors for sample-level statistics
52
+## The signal to noise ratio (SNR)
53
+cnSet$SNR[1:5]
54
+## the skew
55
+cnSet$SKW[1:5]
56
+## the gender (gender is imputed unless specified in the call to crlmm)
57
+table(cnSet$gender)  ## 1=male, 2=female
58
+## --------------------------------------------------
59
+
60
+
61
+## --------------------------------------------------
62
+## accessors for linear model parameters estimated from
63
+## the linear model for copy number
64
+## (note that the parameters have dimension R x C, where
65
+##  R corresponds to the number of features and C corresponds
66
+##  to the number of batches)
67
+## --------------------------------------------------
68
+## estimate of background
69
+dim(nu(cnSet, "A"))
70
+## background for the A allele in the 2 batches for the
71
+## first 5 markers
72
+nu(cnSet, "A")[1:5, ]
73
+## background for the B allele in the 2 batches for the
74
+## first 5 markers
75
+nu(cnSet, "B")[1:5, ]
76
+## the slope
77
+phi(cnSet, "A")[1:5, ]
78
+## the variance for CN > 0 (log2-scale)
79
+sigma2(cnSet, "A")[1:5, ]
80
+sigma2(cnSet, "B")[1:5, ]
81
+## background variance (log2-scale)
82
+tau2(cnSet, "A")[1:5, ]
83
+tau2(cnSet, "B")[1:5, ]
84
+## correlation within genotype cluster AA
85
+corr(cnSet, "AA")[1:5, ]
86
+## correlation within genotype cluster AB
87
+corr(cnSet, "AB")[1:5, ]
88
+## correlation within genotype cluster BB
89
+corr(cnSet, "BB")[1:5, ]
90
+## --------------------------------------------------
91
+
92
+## --------------------------------------------------
93
+## calculating allele-specific copy number
94
+## --------------------------------------------------
95
+## copy number for allele A, first 5 markers, first 2 samples
96
+(ca <- CA(cnSet, i=1:5, j=1:2))
97
+## copy number for allele B, first 5 markers, first 2 samples
98
+(cb <- CB(cnSet, i=1:5, j=1:2))
99
+## total copy number for first 5 markers, first 2 samples
100
+(cn1 <- ca+cb)
101
+
102
+## total copy number at first 5 nonpolymorphic loci
103
+index <- which(!isSnp(cnSet))[1:5]
104
+cn2 <- CA(cnSet, i=index, j=1:2)
105
+## note, cb is NA at nonpolymorphic loci
106
+(cb <- CB(cnSet, i=index, j=1:2))
107
+## note, ca+cb will give NAs at nonpolymorphic loci
108
+CA(cnSet, i=index, j=1:2) + cb
109
+## A shortcut for total copy number
110
+cn3 <- totalCopyNumber(cnSet, i=1:5, j=1:2)
111
+all.equal(cn3, cn1)
112
+cn4 <- totalCopyNumber(cnSet, i=index, j=1:2)
113
+all.equal(cn4, cn2)
114
+
115
+## markers 1-5, all samples
116
+cn5 <- totalCopyNumber(cnSet, i=1:5)
117
+## all markers, samples 1-5
118
+cn6 <- totalCopyNumber(cnSet, j=1:5)
119
+
120
+## NOTE: subsetting the object before extracting copy number
121
+##       can be very inefficient when the data set is very large,
122
+##       particularly if using ff objects.  IN particular, subsetting
123
+##       the CNSet object will subset all of the assay data elements
124
+##       and all of the elements in the LinearModelParameter slot
125
+\dontrun{
126
+	cnsubset <- cnSet[1:5, ]
127
+}
33 128
 }
34 129
 \keyword{datasets}