Browse code

Use original version of bivariateCenter.

Defining the bivariate center to acknowledge background through correlation did not help the signal to noise ratio in the trisomy data. This code is now in .bivariateCenter.

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

Rob Scharp authored on 01/10/2011 04:49:31
Showing2 changed files

... ...
@@ -1,7 +1,7 @@
1 1
 Package: crlmm
2 2
 Type: Package
3 3
 Title: Genotype Calling (CRLMM) and Copy Number Analysis tool for Affymetrix SNP 5.0 and 6.0 and Illumina arrays.
4
-Version: 1.11.32
4
+Version: 1.11.33
5 5
 Date: 2010-12-10
6 6
 Author: Benilton S Carvalho <Benilton.Carvalho@cancer.org.uk>, Robert Scharpf <rscharpf@jhsph.edu>, Matt Ritchie <mritchie@wehi.edu.au>, Ingo Ruczinski <iruczins@jhsph.edu>, Rafael A Irizarry
7 7
 Maintainer: Benilton S Carvalho <Benilton.Carvalho@cancer.org.uk>, Robert Scharpf <rscharpf@jhsph.edu>, Matt Ritchie <mritchie@wehi.EDU.AU>
... ...
@@ -384,6 +384,36 @@ setMethod("calculatePosteriorMean", signature(object="CNSet"),
384 384
 		  return(pm)
385 385
 	  })
386 386
 
387
+		  .bivariateCenter <- function(nu, phi){
388
+			  ##  lexical scope for mus, CA, CB
389
+			  if(CA <= 2 & CB <= 2 & (CA+CB) < 4){
390
+				  mus[,1, ] <- log2(nu[, 1, ] + CA *
391
+						    phi[, 1, ])
392
+				  mus[,2, ] <- log2(nu[, 2, ] + CB *
393
+						    phi[, 2, ])
394
+			  } else { ## CA > 2
395
+				  if(CA > 2){
396
+					  theta <- pi/4*Sigma[,2,]
397
+					  shiftA <- CA/4*phi[, 1, ] * cos(theta)
398
+					  shiftB <- CA/4*phi[, 1, ] * sin(theta)
399
+					  mus[, 1, ] <- log2(nu[, 1, ] + 2 * phi[,1,]+shiftA)
400
+					  mus[, 2, ] <- log2(nu[, 2, ] + CB *phi[,2,]+shiftB)
401
+				  }
402
+				  if(CB > 2){
403
+					  ## CB > 2
404
+					  theta <- pi/2-pi/4*Sigma[,2,]
405
+					  shiftA <- CB/4*phi[, 2, ] * cos(theta)
406
+					  shiftB <- CB/4*phi[, 2, ] * sin(theta)
407
+					  mus[, 1, ] <- log2(nu[, 1, ] + CA*phi[,1,]+shiftA)
408
+					  mus[, 2, ] <- log2(nu[, 2, ]+ 2*phi[,2,]+shiftB)
409
+				  }
410
+				  if(CA == 2 & CB == 2){
411
+					  mus[, 1, ] <- log2(nu[, 1, ] + 1/2*CA*phi[,1,])
412
+					  mus[, 2, ] <- log2(nu[, 2, ]+ 1/2*CB*phi[,2,])
413
+				  }
414
+			  }
415
+			  mus
416
+		  }
387 417
 
388 418
 ## for a given copy number, return a named list of bivariate normal prediction regions
389 419
 ##   - elements of list are named by genotype
... ...
@@ -428,9 +458,11 @@ setMethod("predictionRegion", signature(object="CNSet", copyNumber="integer"),
428 458
 				 bnms))
429 459
 		  bivariateCenter <- function(nu, phi){
430 460
 			  ##  lexical scope for mus, CA, CB
431
-			  mus[,1, ] <- log2(nu[, 1, ] + CA * phi[, 1, ])
432
-			  mus[,2, ] <- log2(nu[, 2, ] + CB * phi[, 2, ])
433
-			  mus
461
+			  mus[,1, ] <- log2(nu[, 1, ] + CA *
462
+					    phi[, 1, ])
463
+			  mus[,2, ] <- log2(nu[, 2, ] + CB *
464
+					    phi[, 2, ])
465
+			  return(mus)
434 466
 		  }
435 467
 		  np.index <- which(!isSnp(object))
436 468
 		  for(i in seq_along(copyNumber)){
... ...
@@ -452,14 +484,20 @@ setMethod("predictionRegion", signature(object="CNSet", copyNumber="integer"),
452 484
 					  Sigma[, 1, ] <- getVar(object, "tau2A.AA")
453 485
 					  Sigma[np.index, 2, ] <- NA
454 486
 				  }
455
-				  res[[gt]]$mu <- bivariateCenter(nu, phi)
487
+				  res[[gt]]$mu <- bivariateCenter(nu,phi)
488
+				  ## adjust correlation
489
+##				  if(CA == 0 | CB == 0){
490
+##					  Sigma[,2,] <- 0
491
+##				  }
456 492
 				  res[[gt]]$cov <- Sigma
457 493
 			  }
458 494
 		  }
495
+		  res <- as(res, "PredictionRegion")
459 496
 		  return(res)
460 497
 	  })
461 498
 
462 499
 
500
+
463 501
 setMethod("xyplotcrlmm", signature(x="formula", data="CNSet", predictRegion="list"),
464 502
 	  function(x, data, predictRegion, ...){
465 503
 		  fns <- featureNames(data)
... ...
@@ -479,7 +517,7 @@ setMethod("xyplotcrlmm", signature(x="formula", data="CNSet", predictRegion="lis
479 517
 			  x$cov <- x$cov[, , batch.index, drop=FALSE]
480 518
 			  return(x)
481 519
 		  }, bns=bns)
482
-		  ##df <- as.data.frame(data)
520
+		  ##predictRegion is an argument of ABpanel
483 521
 		  xyplot(x, df, predictRegion=predictRegion, ...)
484 522
 	  })
485 523
 setMethod("xyplot", signature(x="formula", data="CNSet"),