Browse code

excplicit coercion of A() and B() to matrices in C1, C2, and C3 functions (affect totalCopynumber)

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

Rob Scharp authored on 16/02/2011 15:59:48
Showing2 changed files

... ...
@@ -390,20 +390,6 @@ shrinkGenotypeSummaries <- function(strata, index.list, object, MIN.OBS, MIN.SAM
390 390
 	shrink.corrAB <- corrAB <- as.matrix(corrAB(object)[marker.index, ])
391 391
 	shrink.corrBB <- corrBB <- as.matrix(corrBB(object)[marker.index, ])
392 392
 	flags <- as.matrix(flags(object)[marker.index, ])
393
-	if(length(batches) >= 3){
394
-		if(verbose) message("Imputing centers for unobserved genotypes from other batches")
395
-		res <- imputeAcrossBatch(N.AA, N.AB, N.BB,
396
-					 medianA.AA, medianA.AB, medianA.BB,
397
-					 medianB.AA, medianB.AB, medianB.BB)
398
-		medianA.AA <- res[["medianA.AA"]]
399
-		medianA.AB <- res[["medianA.AB"]]
400
-		medianA.BB <- res[["medianA.BB"]]
401
-		medianB.AA <- res[["medianB.AA"]]
402
-		medianB.AB <- res[["medianB.AB"]]
403
-		medianB.BB <- res[["medianB.BB"]]
404
-		updated <- res[["updated"]]
405
-	}
406
-
407 393
 	for(k in seq(along=batches)){
408 394
 		sample.index <- batches[[k]]
409 395
 		this.batch <- unique(as.character(batch(object)[sample.index]))
... ...
@@ -473,6 +459,19 @@ shrinkGenotypeSummaries <- function(strata, index.list, object, MIN.OBS, MIN.SAM
473 459
 		negB <- rowSums(medianB[[k]] < 0) > 0
474 460
 		flags[, k] <- as.integer(rowSums(NN == 0) > 0 | negA | negB)
475 461
 	}
462
+	if(length(batches) >= 3){
463
+		if(verbose) message("Imputing centers for unobserved genotypes from other batches")
464
+		res <- imputeAcrossBatch(N.AA, N.AB, N.BB,
465
+					 medianA.AA, medianA.AB, medianA.BB,
466
+					 medianB.AA, medianB.AB, medianB.BB)
467
+		medianA.AA <- res[["medianA.AA"]]
468
+		medianA.AB <- res[["medianA.AB"]]
469
+		medianA.BB <- res[["medianA.BB"]]
470
+		medianB.AA <- res[["medianB.AA"]]
471
+		medianB.AB <- res[["medianB.AB"]]
472
+		medianB.BB <- res[["medianB.BB"]]
473
+		updated <- res[["updated"]]
474
+	}
476 475
 	flags(object)[marker.index, ] <- flags
477 476
 	medianA.AA(object)[marker.index, ] <- do.call("cbind", lapply(medianA, function(x) x[, 1]))
478 477
 	medianA.AB(object)[marker.index, ] <- do.call("cbind", lapply(medianA, function(x) x[, 2]))
... ...
@@ -318,7 +318,7 @@ C1 <- function(object, marker.index, batch.index, sample.index){
318 318
 		jj <- sample.index[as.character(batch(object))[sample.index] == batchNames(object)[l]]
319 319
 		bg <- nuA(object)[marker.index, l]
320 320
 		slope <- phiA(object)[marker.index, l]
321
-		I <- A(object)[marker.index, jj]
321
+		I <- as.matrix(A(object)[marker.index, jj])
322 322
 		acn[, match(jj, sample.index)] <- 1/slope * (I - bg)
323 323
 	}
324 324
 ##	if(length(acn) > 1){
... ...
@@ -338,8 +338,8 @@ C2 <- function(object, marker.index, batch.index, sample.index, NP.X=FALSE){
338 338
 		bg <- nuB(object)[marker.index, l]
339 339
 		slope <- phiB(object)[marker.index, l]
340 340
 		if(!NP.X){
341
-			I <- B(object)[marker.index, jj]
342
-		} else I <- A(object)[marker.index, jj]
341
+			I <- as.matrix(B(object)[marker.index, jj])
342
+		} else I <- as.matrix(A(object)[marker.index, jj])
343 343
 		acn[, match(jj, sample.index)] <- 1/slope * (I - bg)
344 344
 	}
345 345
 ##	if(length(acn) > 1){
... ...
@@ -363,10 +363,8 @@ C3 <- function(object, allele, marker.index, batch.index, sample.index){
363 363
 		phiB <- phiB(object)[marker.index, l]
364 364
 		nuA <- nuA(object)[marker.index, l]
365 365
 		nuB <- nuB(object)[marker.index, l]
366
-		IA <- A(object)[marker.index, jj]
367
-		IB <- B(object)[marker.index, jj]
368
-		## I = nu + acn * phi
369
-		## acn = 1/phi* (I-nu)
366
+		IA <- as.matrix(A(object)[marker.index, jj])
367
+		IB <- as.matrix(B(object)[marker.index, jj])
370 368
 		phistar <- phiB2/phiA
371 369
 		tmp <- (IB - nuB - phistar*IA + phistar*nuA)/phiB
372 370
 		CB <- tmp/(1-phistar*phiA2/phiB)