Browse code

posteriorMean.snp replicated scale.sd to length 2 if it was provided as a scalar

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

Rob Scharp authored on 01/10/2011 04:45:16
Showing1 changed files

... ...
@@ -2475,7 +2475,14 @@ calculatePosteriorMean <- function(object, type=c("SNP", "NP", "X.SNP", "X.NP"),
2475 2475
 
2476 2476
 	#}
2477 2477
 	} else stop("type not available")
2478
-	if(length(emit)==1) emit <- emit[[1]] else stop("need to rbind elements of emit list?")
2478
+	if(length(emit)==1) emit <- emit[[1]] else {
2479
+		tmp <- array(NA, dim=c(length(unlist(marker.list)), ncol(object), length(prior.prob)))
2480
+		for(i in seq_along(marker.list)){
2481
+			marker.index <- marker.list[[i]]
2482
+			tmp[marker.index, , ] <- emit[[i]]
2483
+		}
2484
+		emit <- tmp
2485
+	}#stop("need to rbind elements of emit list?")
2479 2486
 	##tmp <- do.call("rbind", emit)
2480 2487
 	match.index <- match(rownames(emit), featureNames(object))
2481 2488
 	S <- length(prior.prob)
... ...
@@ -2485,14 +2492,7 @@ calculatePosteriorMean <- function(object, type=c("SNP", "NP", "X.SNP", "X.NP"),
2485 2492
 	return(object)
2486 2493
 }
2487 2494
 
2488
-posteriorMean.snp <- function(stratum, object, index.list, CN,
2489
-			      prior.prob=c(1/7, 1/7, 3/7, 1/7, 1/7), is.lds=TRUE, verbose=TRUE,
2490
-			      scale.sd=1){
2491
-	if(length(scale.sd) == 1) rep(scale.sd,2)
2492
-	if(verbose) message("Probe stratum ", stratum, " of ", length(index.list))
2493
-	index <- index.list[[stratum]]
2494
-	test <- tryCatch(open(A(object)), error=function(e) NULL)
2495
-	if(!is.null(test)){
2495
+.open <- function(object){
2496 2496
 		open(B(object))
2497 2497
 		open(tau2A.AA(object))
2498 2498
 		open(tau2B.BB(object))
... ...
@@ -2505,7 +2505,47 @@ posteriorMean.snp <- function(stratum, object, index.list, CN,
2505 2505
 		open(nuB(object))
2506 2506
 		open(phiA(object))
2507 2507
 		open(phiB(object))
2508
-	}
2508
+}
2509
+.close <- function(object){
2510
+	close(A(object))
2511
+	close(B(object))
2512
+	close(tau2A.AA(object))
2513
+	close(tau2B.BB(object))
2514
+	close(tau2A.BB(object))
2515
+	close(tau2B.AA(object))
2516
+	close(corrAA(object))
2517
+	close(corrAB(object))
2518
+	close(corrBB(object))
2519
+	close(nuA(object))
2520
+	close(nuB(object))
2521
+	close(phiA(object))
2522
+	close(phiB(object))
2523
+}
2524
+
2525
+sum.mymatrix <- function(...){
2526
+	x <- list(...)
2527
+	return(x[[1]] + do.call(sum, x[-1]))
2528
+}
2529
+numberGenotypes <- function(CT){
2530
+	stopifnot(length(CT)==1)
2531
+	copynumber <- paste("cn", CT, sep="")
2532
+	switch(copynumber,
2533
+		cn0=1,
2534
+		cn1=2,
2535
+		cn2=3,
2536
+		cn3=4,
2537
+		cn4=4,
2538
+		cn5=6, NULL)
2539
+}
2540
+
2541
+posteriorMean.snp <- function(stratum, object, index.list, CN,
2542
+			      prior.prob=c(1/7, 1/7, 3/7, 1/7, 1/7), is.lds=TRUE, verbose=TRUE,
2543
+			      scale.sd=1){
2544
+	if(length(scale.sd) == 1) rep(scale.sd,2)
2545
+	if(verbose) message("Probe stratum ", stratum, " of ", length(index.list))
2546
+	index <- index.list[[stratum]]
2547
+	test <- tryCatch(open(A(object)), error=function(e) NULL)
2548
+	if(!is.null(test)) .open(object)
2509 2549
 	a <- log2(as.matrix(A(object)[index, ]))
2510 2550
 	b <- log2(as.matrix(B(object)[index, ]))
2511 2551
 	NN <- Ns(object, i=index)[, , 1]
... ...
@@ -2520,39 +2560,10 @@ posteriorMean.snp <- function(stratum, object, index.list, CN,
2520 2560
 	phiA <- as.matrix(phiA(object)[index, ])
2521 2561
 	nuB <- as.matrix(nuB(object)[index, ])
2522 2562
 	phiB <- as.matrix(phiB(object)[index, ])
2523
-	if(!is.null(test)){
2524
-		close(A(object))
2525
-		close(B(object))
2526
-		close(tau2A.AA(object))
2527
-		close(tau2B.BB(object))
2528
-		close(tau2A.BB(object))
2529
-		close(tau2B.AA(object))
2530
-		close(corrAA(object))
2531
-		close(corrAB(object))
2532
-		close(corrBB(object))
2533
-		close(nuA(object))
2534
-		close(nuB(object))
2535
-		close(phiA(object))
2536
-		close(phiB(object))
2537
-	}
2563
+	if(!is.null(test)) .close(object)
2538 2564
 	S <- length(prior.prob)
2539 2565
 	emit <- array(NA, dim=c(nrow(a), ncol(a), S))##SNPs x sample x 'truth'
2540 2566
 	sample.index <- split(1:ncol(object), batch(object))
2541
-	sum.mymatrix <- function(...){
2542
-		x <- list(...)
2543
-		return(x[[1]] + do.call(sum, x[-1]))
2544
-	}
2545
-	numberGenotypes <- function(CT){
2546
-		stopifnot(length(CT)==1)
2547
-		copynumber <- paste("cn", CT, sep="")
2548
-		switch(copynumber,
2549
-		       cn0=1,
2550
-		       cn1=2,
2551
-		       cn2=3,
2552
-		       cn3=4,
2553
-		       cn4=4,
2554
-		       cn5=6, NULL)
2555
-	}
2556 2567
 	##emit <- vector("list", length(sample.index))
2557 2568
 	for(j in seq_along(sample.index)){
2558 2569
 		cat("batch ", j, "\n")
... ...
@@ -2577,10 +2588,11 @@ posteriorMean.snp <- function(stratum, object, index.list, CN,
2577 2588
 				CB <- CT-CA
2578 2589
 				A.scale <- sqrt(tau2A*(CA==0) + sig2A*(CA > 0))
2579 2590
 				B.scale <- sqrt(tau2B*(CB==0) + sig2B*(CB > 0))
2580
-				if(CA == 0 | CB == 0){
2591
+				if(CA == 0 | CB == 0){## possibly homozygous BB and AA, respectively
2581 2592
 					A.scale <- A.scale*scale.sd[1]
2582 2593
 					B.scale <- B.scale*scale.sd[1]
2583 2594
 				} else { ## one or both greater than zero
2595
+					if(length(scale.sd) == 1) scale.sd <- rep(scale.sd, 2)
2584 2596
 					A.scale <- A.scale*scale.sd[2]
2585 2597
 					B.scale <- B.scale*scale.sd[2]
2586 2598
 				}