Browse code

Commented out computation of standard errors for coefficients. Replaced crossprod(solve...) with call to the Fortran library dqrls

the copynumber vignette has some experimental code that needs to be removed

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

Rob Scharp authored on 21/08/2010 02:47:35
Showing4 changed files

... ...
@@ -73,7 +73,8 @@ export(crlmm,
73 73
        crlmmCopynumber2, crlmmCopynumberLD)
74 74
 export(constructIlluminaCNSet)
75 75
 export(linesCNSetLM)
76
-export(linearModel.noweights, computeCN, fit.lm1, fit.lm2, fit.lm3, fit.lm4, construct)
76
+export(linearModel.noweights, computeCN, fit.lm1, fit.lm2, fit.lm3, fit.lm4, construct,
77
+       dqrlsWrapper, nuphiAllele)
77 78
 
78 79
 
79 80
 
... ...
@@ -8,6 +8,7 @@ setMethod("initialize", "CNSetLM", function(.Object, lM=new("list"), ...){
8 8
 	.Object@lM <- lM
9 9
 	.Object <- callNextMethod(.Object, ...)
10 10
 })
11
+
11 12
 setValidity("CNSetLM",
12 13
 	    function(object){
13 14
 		    if(!"batch" %in% varLabels(protocolData(object)))
... ...
@@ -20,3 +21,6 @@ setValidity("CNSetLM",
20 21
 			    return("'isSnp' not defined in featureData")
21 22
 		    return(TRUE)
22 23
 	    })
24
+
25
+##setClass("CNSetTest", contains="SnpSuperSet")
26
+
... ...
@@ -367,14 +367,18 @@ corByGenotype <- function(A, B, G, Ns, which.cluster=c(1,2,3)[1], DF.PRIOR){
367 367
 	cors
368 368
 }
369 369
 
370
-generateX <- function(w, X) as.numeric(diag(w) %*% X)
371
-generateIXTX <- function(x, nrow=3) {
372
-	X <- matrix(x, nrow=nrow)
373
-	XTX <- crossprod(X)
374
-	solve(XTX)
370
+dqrlsWrapper <- function(x, y, wts, tol=1e-7){
371
+	n <- NROW(y)
372
+	p <- ncol(x)
373
+	ny <- NCOL(y)
374
+	.Fortran("dqrls", qr=x*wts, n=n, p=p, y=y * wts, ny=ny,
375
+		 tol=as.double(tol), coefficients=mat.or.vec(p, ny),
376
+		 residuals=y, effects=mat.or.vec(n, ny),
377
+		 rank=integer(1L), pivot=1L:p, qraux=double(p),
378
+		 work=double(2 * p), PACKAGE="base")[["coefficients"]]
375 379
 }
376 380
 
377
-nuphiAllele2 <- function(allele, Ystar, W, Ns){
381
+fit.wls <- function(allele, Ystar, W, Ns, autosome=TRUE){
378 382
 	complete <- rowSums(is.na(W)) == 0 
379 383
 	if(any(!is.finite(W))){## | any(!is.finite(V))){
380 384
 		i <- which(rowSums(!is.finite(W)) > 0)
... ...
@@ -382,21 +386,23 @@ nuphiAllele2 <- function(allele, Ystar, W, Ns){
382 386
 	}
383 387
 	NOHET <- mean(Ns[, 2], na.rm=TRUE) < 0.05
384 388
 	if(missing(allele)) stop("must specify allele")
385
-	if(allele == "A") X <- cbind(1, 2:0) else X <- cbind(1, 0:2)
389
+	if(autosome){
390
+		if(allele == "A") X <- cbind(1, 2:0)
391
+		if(allele == "B") X <- cbind(1, 0:2)
392
+		betahat <- matrix(NA, 2, nrow(Ystar))
393
+	} else {
394
+		if(allele == "A") X <- cbind(1, c(1, 0, 2, 1, 0), c(0, 1, 0, 1, 2))
395
+		if(allele == "B") X <- cbind(1, c(0, 1, 0, 1, 2), c(1, 0, 2, 1, 0))
396
+		betahat <- matrix(NA, 3, nrow(Ystar))
397
+	}
386 398
 	if(NOHET) X <- X[-2, ] ##more than 1 X chromosome, but all homozygous		
387 399
 	##How to quickly generate Xstar, Xstar = diag(W) %*% X
388
-	Xstar <- apply(W, 1, generateX, X)
389
-	IXTX <- apply(Xstar, 2, generateIXTX, nrow=nrow(X))
390
-	betahat <- matrix(NA, 2, nrow(Ystar))
391
-	ses <- matrix(NA, 2, nrow(Ystar))
400
+	##Xstar <- apply(W, 1, generateX, X)
401
+	ww <- rep(1, ncol(Ystar))
392 402
 	for(i in 1:nrow(Ystar)){
393
-		betahat[, i] <- crossprod(matrix(IXTX[, i], ncol(X), ncol(X)), crossprod(matrix(Xstar[, i], nrow=nrow(X)), Ystar[i, ]))
403
+		betahat[, i] <- dqrlsWrapper(W[i, ] * X, Ystar[i, ], ww)
394 404
 		##ssr <- sum((Ystar[i, ] - matrix(Xstar[, i], nrow(X), ncol(X)) %*% matrix(betahat[, i], ncol(X), 1))^2)
395
-		##ses[, i] <- sqrt(diag(matrix(IXTX[, i], ncol(X), ncol(X)) * ssr))
396 405
 	}
397
-##	nu <- betahat[1, ]
398
-##	phi <- betahat[2, ]
399
-##	return(list(nu, phi))
400 406
 	return(betahat)
401 407
 }
402 408
 
... ...
@@ -426,35 +432,35 @@ linearModel.noweights <- function(allele, Ystar, W, Ns){
426 432
 	return(list(nu, phi))
427 433
 }
428 434
 
429
-nuphiAlleleX <- function(allele, Ystar, W, Ns, chrX=FALSE){
430
-	complete <- rowSums(is.na(W)) == 0 
431
-	if(any(!is.finite(W))){## | any(!is.finite(V))){
432
-		i <- which(rowSums(!is.finite(W)) > 0)
433
-		stop("Possible zeros in the within-genotype estimates of the spread (vA, vB). ")
434
-	}
435
-	##NOHET <- mean(Ns[, 2], na.rm=TRUE) < 0.05
436
-	if(missing(allele)) stop("must specify allele")
437
-	if(allele == "A") X <- cbind(1, c(1, 0, 2, 1, 0), c(0, 1, 0, 1, 2))
438
-	if(allele == "B") X <- cbind(1, c(0, 1, 0, 1, 2), c(1, 0, 2, 1, 0))	
439
-	##if(NOHET) X <- X[-2, ] ##more than 1 X chromosome, but all homozygous		
440
-	##How to quickly generate Xstar, Xstar = diag(W) %*% X
441
-	Xstar <- apply(W, 1, generateX, X)
442
-	IXTX <- apply(Xstar, 2, generateIXTX, nrow=nrow(X))
443
-
444
-	betahat <- matrix(NA, 3, nrow(Ystar))
445
-	ses <- matrix(NA, 3, nrow(Ystar))			
446
-	##betahat <- matrix(NA, 2, nrow(Ystar))
447
-	##ses <- matrix(NA, 2, nrow(Ystar))
448
-	for(i in 1:nrow(Ystar)){
449
-		betahat[, i] <- crossprod(matrix(IXTX[, i], ncol(X), ncol(X)), crossprod(matrix(Xstar[, i], nrow=nrow(X)), Ystar[i, ]))
450
-		ssr <- sum((Ystar[i, ] - matrix(Xstar[, i], nrow(X), ncol(X)) %*% matrix(betahat[, i], ncol(X), 1))^2)
451
-		ses[, i] <- sqrt(diag(matrix(IXTX[, i], ncol(X), ncol(X)) * ssr))
452
-	}
453
-	nu <- betahat[1, ]
454
-	phi <- betahat[2, ]
455
-	phi2 <- betahat[3, ]
456
-	return(list(nu, phi, phi2))
457
-}
435
+##nuphiAlleleX <- function(allele, Ystar, W, Ns, chrX=FALSE){
436
+##	complete <- rowSums(is.na(W)) == 0 
437
+##	if(any(!is.finite(W))){## | any(!is.finite(V))){
438
+##		i <- which(rowSums(!is.finite(W)) > 0)
439
+##		stop("Possible zeros in the within-genotype estimates of the spread (vA, vB). ")
440
+##	}
441
+##	##NOHET <- mean(Ns[, 2], na.rm=TRUE) < 0.05
442
+##	if(missing(allele)) stop("must specify allele")
443
+##	if(allele == "A") X <- cbind(1, c(1, 0, 2, 1, 0), c(0, 1, 0, 1, 2))
444
+##	if(allele == "B") X <- cbind(1, c(0, 1, 0, 1, 2), c(1, 0, 2, 1, 0))	
445
+##	##if(NOHET) X <- X[-2, ] ##more than 1 X chromosome, but all homozygous		
446
+##	##How to quickly generate Xstar, Xstar = diag(W) %*% X
447
+##	Xstar <- apply(W, 1, generateX, X)
448
+##	IXTX <- apply(Xstar, 2, generateIXTX, nrow=nrow(X))
449
+##
450
+##	betahat <- matrix(NA, 3, nrow(Ystar))
451
+##	ses <- matrix(NA, 3, nrow(Ystar))			
452
+##	##betahat <- matrix(NA, 2, nrow(Ystar))
453
+##	##ses <- matrix(NA, 2, nrow(Ystar))
454
+##	for(i in 1:nrow(Ystar)){
455
+##		betahat[, i] <- crossprod(matrix(IXTX[, i], ncol(X), ncol(X)), crossprod(matrix(Xstar[, i], nrow=nrow(X)), Ystar[i, ]))
456
+##		ssr <- sum((Ystar[i, ] - matrix(Xstar[, i], nrow(X), ncol(X)) %*% matrix(betahat[, i], ncol(X), 1))^2)
457
+##		ses[, i] <- sqrt(diag(matrix(IXTX[, i], ncol(X), ncol(X)) * ssr))
458
+##	}
459
+##	nu <- betahat[1, ]
460
+##	phi <- betahat[2, ]
461
+##	phi2 <- betahat[3, ]
462
+##	return(list(nu, phi, phi2))
463
+##}
458 464
 
459 465
 
460 466
 nuphiAllele <- function(object, allele, Ystar, W, tmp.objects, cnOptions){
... ...
@@ -470,18 +476,17 @@ nuphiAllele <- function(object, allele, Ystar, W, tmp.objects, cnOptions){
470 476
 	}	
471 477
 	Ns <- tmp.objects[["Ns"]]	
472 478
 	Ns <- Ns[I, ]
473
-	
474 479
 	CHR <- unique(chromosome(object))
475 480
 	##batch <- unique(object$batch)
476 481
 	batch <- unique(batch(object))
477 482
 	nuA <- getParam(object, "nuA", batch)
478 483
 	nuB <- getParam(object, "nuB", batch)
479
-	nuA.se <- getParam(object, "nuA.se", batch)
480
-	nuB.se <- getParam(object, "nuB.se", batch)
484
+##	nuA.se <- getParam(object, "nuA.se", batch)
485
+##	nuB.se <- getParam(object, "nuB.se", batch)
481 486
 	phiA <- getParam(object, "phiA", batch)
482 487
 	phiB <- getParam(object, "phiB", batch)
483
-	phiA.se <- getParam(object, "phiA.se", batch)
484
-	phiB.se <- getParam(object, "phiB.se", batch)	
488
+##	phiA.se <- getParam(object, "phiA.se", batch)
489
+##	phiB.se <- getParam(object, "phiB.se", batch)	
485 490
 	if(CHR==23){
486 491
 		phiAX <- getParam(object, "phiAX", batch)
487 492
 		phiBX <- getParam(object, "phiBX", batch)		
... ...
@@ -494,38 +499,36 @@ nuphiAllele <- function(object, allele, Ystar, W, tmp.objects, cnOptions){
494 499
 			if(allele == "A") X <- cbind(1, c(1, 0, 2, 1, 0), c(0, 1, 0, 1, 2))
495 500
 			if(allele == "B") X <- cbind(1, c(0, 1, 0, 1, 2), c(1, 0, 2, 1, 0))
496 501
 		} else{
497
-			if(allele == "A") X <- cbind(1, c(1, 0, 2, 0), c(0, 1, 0, 2))
502
+		if(allele == "A") X <- cbind(1, c(1, 0, 2, 0), c(0, 1, 0, 2))
498 503
 			if(allele == "B") X <- cbind(1, c(0, 1, 0, 2), c(1, 0, 2, 0))
499 504
 		}
505
+		betahat <- matrix(NA, 3, nrow(Ystar))
500 506
 	} else {##autosome
501 507
 		if(allele == "A") X <- cbind(1, 2:0) else X <- cbind(1, 0:2)
502
-		if(NOHET) X <- X[-2, ] ##more than 1 X chromosome, but all homozygous		
508
+		if(NOHET) X <- X[-2, ] ##more than 1 X chromosome, but all homozygous
509
+		betahat <- matrix(NA, 2, nrow(Ystar))
503 510
 	}
504 511
 	##How to quickly generate Xstar, Xstar = diag(W) %*% X
505
-	Xstar <- apply(W, 1, generateX, X)
506
-	IXTX <- apply(Xstar, 2, generateIXTX, nrow=nrow(X))
512
+##	Xstar <- apply(W, 1, generateX, X)
513
+##	IXTX <- apply(Xstar, 2, generateIXTX, nrow=nrow(X))
507 514
 	##as.numeric(diag(W[1, ]) %*% X)
508
-	if(CHR == 23){
509
-		betahat <- matrix(NA, 3, nrow(Ystar))
510
-		ses <- matrix(NA, 3, nrow(Ystar))		
511
-	} else{
512
-		betahat <- matrix(NA, 2, nrow(Ystar))
513
-		ses <- matrix(NA, 2, nrow(Ystar))
514
-	}
515
-	for(i in 1:nrow(Ystar)){
516
-		betahat[, i] <- crossprod(matrix(IXTX[, i], ncol(X), ncol(X)), crossprod(matrix(Xstar[, i], nrow=nrow(X)), Ystar[i, ]))
517
-		ssr <- sum((Ystar[i, ] - matrix(Xstar[, i], nrow(X), ncol(X)) %*% matrix(betahat[, i], ncol(X), 1))^2)
518
-		ses[, i] <- sqrt(diag(matrix(IXTX[, i], ncol(X), ncol(X)) * ssr))
515
+	ww <- rep(1, ncol(Ystar))
516
+	II <- which(rowSums(is.nan(Ystar)) == 0)
517
+	for(i in II){
518
+		betahat[, i] <- dqrlsWrapper(W[i, ] * X, Ystar[i, ], ww)
519
+		##betahat[, i] <- crossprod(matrix(IXTX[, i], ncol(X), ncol(X)), crossprod(matrix(Xstar[, i], nrow=nrow(X)), Ystar[i, ]))
520
+		##ssr <- sum((Ystar[i, ] - matrix(Xstar[, i], nrow(X), ncol(X)) %*% matrix(betahat[, i], ncol(X), 1))^2)
521
+		##ses[, i] <- sqrt(diag(matrix(IXTX[, i], ncol(X), ncol(X)) * ssr))
519 522
 	}
520 523
 	if(allele == "A"){
521 524
 		nuA[complete] <- betahat[1, ]
522 525
 		phiA[complete] <- betahat[2, ]
523
-		nuA.se[complete] <- ses[1, ]
524
-		phiA.se[complete] <- ses[2, ]
526
+##		nuA.se[complete] <- ses[1, ]
527
+##		phiA.se[complete] <- ses[2, ]
525 528
 		object <- pr(object, "nuA", batch, nuA)
526
-		object <- pr(object, "nuA.se", batch, nuA.se)
529
+##		object <- pr(object, "nuA.se", batch, nuA.se)
527 530
 		object <- pr(object, "phiA", batch, phiA)
528
-		object <- pr(object, "phiA.se", batch, phiA.se)
531
+##		object <- pr(object, "phiA.se", batch, phiA.se)
529 532
 		if(CHR == 23){
530 533
 			phiAX[complete] <- betahat[3, ]
531 534
 			object <- pr(object, "phiAX", batch, phiAX)			
... ...
@@ -534,12 +537,12 @@ nuphiAllele <- function(object, allele, Ystar, W, tmp.objects, cnOptions){
534 537
 	if(allele=="B"){
535 538
 		nuB[complete] <- betahat[1, ]
536 539
 		phiB[complete] <- betahat[2, ]
537
-		nuB.se[complete] <- ses[1, ]
538
-		phiB.se[complete] <- ses[2, ]
540
+##		nuB.se[complete] <- ses[1, ]
541
+##		phiB.se[complete] <- ses[2, ]
539 542
 		object <- pr(object, "nuB", batch, nuB)
540
-		object <- pr(object, "nuB.se", batch, nuB.se)
543
+##		object <- pr(object, "nuB.se", batch, nuB.se)
541 544
 		object <- pr(object, "phiB", batch, phiB)		
542
-		object <- pr(object, "phiB.se", batch, phiB.se)
545
+##		object <- pr(object, "phiB.se", batch, phiB.se)
543 546
 		if(CHR == 23){
544 547
 			phiBX[complete] <- betahat[3, ]
545 548
 			object <- pr(object, "phiBX", batch, phiBX)
... ...
@@ -866,7 +869,9 @@ fit.lm1 <- function(idxBatch,
866 869
 	CP <- as.matrix(snpCallProbability(object)[snps, ])
867 870
 	AA <- as.matrix(A(object)[snps, ])
868 871
 	BB <- as.matrix(B(object)[snps, ])
872
+	zzzz <- 1
869 873
 	for(k in batches){
874
+		zzzz <- zzzz+1
870 875
 		G <- GG[, k]
871 876
 		NORM <- normal.snps[, k]
872 877
 		xx <- CP[, k]
... ...
@@ -957,13 +962,19 @@ fit.lm1 <- function(idxBatch,
957 962
 		YA <- muA*wA
958 963
 		YB <- muB*wB
959 964
 		if(weighted.lm){
960
-			res <- nuphiAllele2(allele="A", Ystar=YA, W=wA, Ns=Ns)
961
-		} else res <- linearModel.noweights(allele="A", Ystar=YA, W=wA, Ns=Ns)
965
+			res <- fit.wls(allele="A", Ystar=YA, W=wA, Ns=Ns)
966
+		} else{
967
+			if(zzzz==1) message("currently, only weighted least squares (wls) is available... fitting wls")
968
+			res <- fit.wls(allele="A", Ystar=YA, W=wA, Ns=Ns)			
969
+		}
962 970
 		nuA[, J] <- res[[1]]
963 971
 		phiA[, J] <- res[[2]]
964
-		if(!weighted.lm){
965
-			res <- nuphiAllele2(allele="B", Ystar=YB, W=wB, Ns=Ns)
966
-		} else res <- linearModel.noweights(allele="B", Ystar=YB, W=wB, Ns=Ns)
972
+		if(weighted.lm){
973
+			res <- fit.wls(allele="B", Ystar=YB, W=wB, Ns=Ns)
974
+		} else {
975
+			if(zzzz==1) message("currently, only weighted least squares (wls) is available... fitting wls")
976
+			res <- fit.wls(allele="B", Ystar=YB, W=wB, Ns=Ns)			
977
+		}
967 978
 		##nuB[, J] <- res[[1]]
968 979
 		nuB[, J] <- res[1, ]
969 980
 		##phiB[, J] <- res[[2]]
... ...
@@ -1137,7 +1148,6 @@ fit.lm2 <- function(idxBatch,
1137 1148
 	cA[cA > 5] <-  5
1138 1149
 	cA <- matrix(as.integer(cA*100), nrow(cA), ncol(cA))
1139 1150
 	CA(object)[snps, ] <- cA
1140
-	##open(lM(object))
1141 1151
 	tmp <- physical(lM(object))$nuA
1142 1152
 	tmp[snps, ] <- nuA.np
1143 1153
 	lM(object)$nuA <- tmp
... ...
@@ -1147,9 +1157,6 @@ fit.lm2 <- function(idxBatch,
1147 1157
 	tmp <- physical(lM(object))$phiA
1148 1158
 	tmp[snps, ] <- phiA.np
1149 1159
 	lM(object)$sig2A <- tmp
1150
-##	lM(object)$sig2A[snps, ] <- sig2A.np
1151
-##	lM(object)$nuA[snps, ] <- nuA.np
1152
-##	lM(object)$phiA[snps, ] <- phiA.np
1153 1160
 	lapply(assayData(object), close)
1154 1161
 	lapply(lM(object), close)
1155 1162
 	TRUE
... ...
@@ -1306,14 +1313,16 @@ fit.lm3 <- function(idxBatch,
1306 1313
 		wB <- sqrt(1/vB2)
1307 1314
 		YA <- muA*wA
1308 1315
 		YB <- muB*wB
1309
-		res <- nuphiAlleleX(allele="A", Ystar=YA, W=wA)
1310
-		nuA[, J] <- res[[1]]
1311
-		phiA[, J] <- res[[2]]
1312
-		phiA2[, J] <- res[[3]]
1313
-		res <- nuphiAlleleX(allele="B", Ystar=YB, W=wB)
1314
-		nuB[, J] <- res[[1]]
1315
-		phiB[, J] <- res[[2]]
1316
-		phiB2[, J] <- res[[3]]
1316
+		##res <- nuphiAlleleX(allele="A", Ystar=YA, W=wA)
1317
+		betas <- fit.wls(allele="A", Ystar=YA, W=wA, autosome=FALSE)
1318
+		nuA[, J] <- betas[1, ]
1319
+		phiA[, J] <- betas[2, ]
1320
+		phiA2[, J] <- betas[3, ]
1321
+		rm(betas)
1322
+		betas <- fit.wls(allele="B", Ystar=YB, W=wB, autosome=FALSE)
1323
+		nuB[, J] <- betas[1, ]
1324
+		phiB[, J] <- betas[2, ]
1325
+		phiB2[, J] <- betas[3, ]
1317 1326
 		if(THR.NU.PHI){
1318 1327
 			nuA[nuA[, J] < MIN.NU, J] <- MIN.NU
1319 1328
 			nuB[nuB[, J] < MIN.NU, J] <- MIN.NU
... ...
@@ -2715,74 +2724,74 @@ thresholdModelParams <- function(object, cnOptions){
2715 2724
 ##}
2716 2725
 
2717 2726
 
2718
-##computeCopynumber.CNSet <- function(object, cnOptions){
2719
-##	##PLATE <- unique(object$batch)
2720
-##	PLATE <- unique(batch(object))
2721
-##	verbose <- cnOptions$verbose
2722
-##	tmp.objects <- instantiateObjects(object, cnOptions)
2723
-##	bias.adj <- cnOptions$bias.adj
2724
-##	if(bias.adj & ncol(object) <= 15){
2725
-##		warning(paste("bias.adj is TRUE, but too few samples to perform this step"))
2726
-##		cnOptions$bias.adj <- bias.adj <- FALSE
2727
-##	}
2728
-##	if(bias.adj){
2729
-##		if(verbose) message("Dropping samples with low posterior prob. of normal copy number (samples dropped is locus-specific)")
2730
-##		tmp.objects <- biasAdjNP(object, cnOptions, tmp.objects)
2731
-##		tmp.objects <- biasAdj(object, cnOptions, tmp.objects)
2732
-##		if(verbose) message("Recomputing location and scale parameters")
2733
-##	}
2734
-##	##update tmp.objects
2735
-##	tmp.objects <- withinGenotypeMoments(object,
2736
-##					     cnOptions=cnOptions,
2737
-##					     tmp.objects=tmp.objects)
2738
-##	object <- locationAndScale(object, cnOptions, tmp.objects)
2739
-##	tmp.objects <- oneBatch(object,
2740
-##				cnOptions=cnOptions,
2741
-##				tmp.objects=tmp.objects)
2742
-##	##coefs calls nuphiAllele.
2743
-##	object <- coefs(object, cnOptions, tmp.objects)
2744
-##	##nuA=getParam(object, "nuA", PLATE)
2745
-##	THR.NU.PHI <- cnOptions$THR.NU.PHI
2746
-##	if(THR.NU.PHI){
2747
-##		verbose <- cnOptions$verbose
2748
-##		##if(verbose) message("Thresholding nu and phi")
2749
-##		object <- thresholdModelParams(object, cnOptions)
2750
-##	}		
2751
-##	##if(verbose) message("\nAllele specific copy number")	
2752
-##	object <- polymorphic(object, cnOptions, tmp.objects)
2753
-##	if(any(!isSnp(object))){ ## there are nonpolymorphic probes
2754
-##		##if(verbose) message("\nCopy number for nonpolymorphic probes...")	
2755
-##		object <- nonpolymorphic(object, cnOptions, tmp.objects)
2756
-##	}
2757
-##	##---------------------------------------------------------------------------
2758
-##	##Note: the replacement method multiples by 100
2759
-####	CA(object)[, batch==PLATE] <- CA(object)
2760
-####	CB(object)[, batch==PLATE] <- CB(object)
2761
-##	##---------------------------------------------------------------------------
2762
-##	##update-the plate-specific parameters for copy number
2763
-##	object <- pr(object, "nuA", PLATE, getParam(object, "nuA", PLATE))
2764
-##	object <- pr(object, "nuA.se", PLATE, getParam(object, "nuA.se", PLATE))
2765
-##	object <- pr(object, "nuB", PLATE, getParam(object, "nuB", PLATE))
2766
-##	object <- pr(object, "nuB.se", PLATE, getParam(object, "nuB.se", PLATE))
2767
-##	object <- pr(object, "phiA", PLATE, getParam(object, "phiA", PLATE))
2768
-##	object <- pr(object, "phiA.se", PLATE, getParam(object, "phiA.se", PLATE))
2769
-##	object <- pr(object, "phiB", PLATE, getParam(object, "phiB", PLATE))
2770
-##	object <- pr(object, "phiB.se", PLATE, getParam(object, "phiB.se", PLATE))
2771
-##	object <- pr(object, "tau2A", PLATE, getParam(object, "tau2A", PLATE))
2772
-##	object <- pr(object, "tau2B", PLATE, getParam(object, "tau2B", PLATE))				
2773
-##	object <- pr(object, "sig2A", PLATE, getParam(object, "sig2A", PLATE))
2774
-##	object <- pr(object, "sig2B", PLATE, getParam(object, "sig2B", PLATE))		
2775
-##	object <- pr(object, "phiAX", PLATE, as.numeric(getParam(object, "phiAX", PLATE)))
2776
-##	object <- pr(object, "phiBX", PLATE, as.numeric(getParam(object, "phiBX", PLATE)))
2777
-##	object <- pr(object, "corr", PLATE, getParam(object, "corr", PLATE))
2778
-##	object <- pr(object, "corrA.BB", PLATE, getParam(object, "corrA.BB", PLATE))
2779
-##	object <- pr(object, "corrB.AA", PLATE, getParam(object, "corrB.AA", PLATE))
2780
-##	##object <- object[order(chromosome(object), position(object)), ]
2781
-##	if(cnOptions[["thresholdCopynumber"]]){
2782
-##		object <- thresholdCopynumber(object)
2783
-##	}
2784
-##	return(object)
2785
-##}
2727
+computeCopynumber.CNSet <- function(object, cnOptions){
2728
+	##PLATE <- unique(object$batch)
2729
+	PLATE <- unique(batch(object))
2730
+	verbose <- cnOptions$verbose
2731
+	tmp.objects <- instantiateObjects(object, cnOptions)
2732
+	bias.adj <- cnOptions$bias.adj
2733
+	if(bias.adj & ncol(object) <= 15){
2734
+		warning(paste("bias.adj is TRUE, but too few samples to perform this step"))
2735
+		cnOptions$bias.adj <- bias.adj <- FALSE
2736
+	}
2737
+	if(bias.adj){
2738
+		if(verbose) message("Dropping samples with low posterior prob. of normal copy number (samples dropped is locus-specific)")
2739
+		tmp.objects <- biasAdjNP(object, cnOptions, tmp.objects)
2740
+		tmp.objects <- biasAdj(object, cnOptions, tmp.objects)
2741
+		if(verbose) message("Recomputing location and scale parameters")
2742
+	}
2743
+	##update tmp.objects
2744
+	tmp.objects <- withinGenotypeMoments(object,
2745
+					     cnOptions=cnOptions,
2746
+					     tmp.objects=tmp.objects)
2747
+	object <- locationAndScale(object, cnOptions, tmp.objects)
2748
+	tmp.objects <- oneBatch(object,
2749
+				cnOptions=cnOptions,
2750
+				tmp.objects=tmp.objects)
2751
+	##coefs calls nuphiAllele.
2752
+	object <- coefs(object, cnOptions, tmp.objects)
2753
+	##nuA=getParam(object, "nuA", PLATE)
2754
+	THR.NU.PHI <- cnOptions$THR.NU.PHI
2755
+	if(THR.NU.PHI){
2756
+		verbose <- cnOptions$verbose
2757
+		##if(verbose) message("Thresholding nu and phi")
2758
+		object <- thresholdModelParams(object, cnOptions)
2759
+	}		
2760
+	##if(verbose) message("\nAllele specific copy number")	
2761
+	object <- polymorphic(object, cnOptions, tmp.objects)
2762
+	if(any(!isSnp(object))){  ##there are nonpolymorphic probes
2763
+		##if(verbose) message("\nCopy number for nonpolymorphic probes...")	
2764
+		object <- nonpolymorphic(object, cnOptions, tmp.objects)
2765
+	}
2766
+	##---------------------------------------------------------------------------
2767
+	##Note: the replacement method multiples by 100
2768
+##	CA(object)[, batch==PLATE] <- CA(object)
2769
+##	CB(object)[, batch==PLATE] <- CB(object)
2770
+	##---------------------------------------------------------------------------
2771
+	##update-the plate-specific parameters for copy number
2772
+	object <- pr(object, "nuA", PLATE, getParam(object, "nuA", PLATE))
2773
+	object <- pr(object, "nuA.se", PLATE, getParam(object, "nuA.se", PLATE))
2774
+	object <- pr(object, "nuB", PLATE, getParam(object, "nuB", PLATE))
2775
+	object <- pr(object, "nuB.se", PLATE, getParam(object, "nuB.se", PLATE))
2776
+	object <- pr(object, "phiA", PLATE, getParam(object, "phiA", PLATE))
2777
+	object <- pr(object, "phiA.se", PLATE, getParam(object, "phiA.se", PLATE))
2778
+	object <- pr(object, "phiB", PLATE, getParam(object, "phiB", PLATE))
2779
+	object <- pr(object, "phiB.se", PLATE, getParam(object, "phiB.se", PLATE))
2780
+	object <- pr(object, "tau2A", PLATE, getParam(object, "tau2A", PLATE))
2781
+	object <- pr(object, "tau2B", PLATE, getParam(object, "tau2B", PLATE))				
2782
+	object <- pr(object, "sig2A", PLATE, getParam(object, "sig2A", PLATE))
2783
+	object <- pr(object, "sig2B", PLATE, getParam(object, "sig2B", PLATE))		
2784
+	object <- pr(object, "phiAX", PLATE, as.numeric(getParam(object, "phiAX", PLATE)))
2785
+	object <- pr(object, "phiBX", PLATE, as.numeric(getParam(object, "phiBX", PLATE)))
2786
+	object <- pr(object, "corr", PLATE, getParam(object, "corr", PLATE))
2787
+	object <- pr(object, "corrA.BB", PLATE, getParam(object, "corrA.BB", PLATE))
2788
+	object <- pr(object, "corrB.AA", PLATE, getParam(object, "corrB.AA", PLATE))
2789
+	##object <- object[order(chromosome(object), position(object)), ]
2790
+	if(cnOptions[["thresholdCopynumber"]]){
2791
+		object <- thresholdCopynumber(object)
2792
+	}
2793
+	return(object)
2794
+}
2786 2795
 
2787 2796
 
2788 2797
 
... ...
@@ -140,6 +140,11 @@ will load \Robject{cnSet.assayData_matrix} from disk if this
140 140
 computation had already been performed as part of the batch job.
141 141
 
142 142
 <<copynumber>>=
143
+stop()
144
+tms1 <- system.time(tmp <- crlmmCopynumber(gtSet.assayData_matrix))
145
+reload_pkg("crlmm")
146
+trace(nuphiAllele, browser)
147
+tms2 <- system.time(tmp <- crlmmCopynumber(gtSet.assayData_matrix))
143 148
 cnSet.assayData_matrix <- checkExists("cnSet.assayData_matrix",
144 149
 				      .path=outdir,
145 150
 				      .FUN=crlmmCopynumber,