Browse code

Update predictionRegion to return arrays in the mu and cov list elements (prediction regions are batch-specific)

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

Rob Scharp authored on 01/10/2011 04:47:24
Showing 4 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.14
4
+Version: 1.11.15
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>
... ...
@@ -102,3 +102,5 @@ setGeneric("calculatePosteriorMean", function(object, posteriorProb, copyNumber=
102 102
 
103 103
 setGeneric("predictionRegion", function(object, copyNumber=0:4)
104 104
 	   standardGeneric("predictionRegion"))
105
+setGeneric("xyplot", useAsDefault=function(x, data, ...) lattice::xyplot(x, data,...))
106
+setGeneric("xyplotcrlmm", function(x, data, predictRegion,...) standardGeneric("xyplotcrlmm"))
... ...
@@ -1,13 +1,4 @@
1
-getNu <- function(object){
2
-	nu.a <- nuA(object)
3
-	nu.b <- nuB(object)
4
-	cbind(nu.a, nu.b)
5
-}
6
-getPhi <- function(object){
7
-	phi.a <- phiA(object)
8
-	phi.b <- phiB(object)
9
-	cbind(phi.a, phi.b)
10
-}
1
+
11 2
 
12 3
 getSigma2 <- function(object){
13 4
 	phi.a <- phiA(object)
... ...
@@ -17,7 +8,7 @@ getSigma2 <- function(object){
17 8
 
18 9
 getCor <- function(object, gt){
19 10
 	if(gt =="NULL") return(rep(0, nrow(object)))
20
-	as.numeric(assayDataElement(batchStatistics(object), paste("corr", gt, sep="")))
11
+	assayDataElement(batchStatistics(object), paste("corr", gt, sep=""))
21 12
 }
22 13
 
23 14
 getTau2 <- function(object, gt){
... ...
@@ -31,7 +22,7 @@ getTau2 <- function(object, gt){
31 22
 
32 23
 getVar <- function(object, nm){
33 24
 	bs <- batchStatistics(object)
34
-	as.numeric(assayDataElement(bs, nm))
25
+	assayDataElement(bs, nm)
35 26
 }
36 27
 
37 28
 setMethod("nuA", signature=signature(object="CNSet"), function(object) nu(object, "A"))
... ...
@@ -345,19 +345,38 @@ setMethod("predictionRegion", signature(object="CNSet", copyNumber="integer"),
345 345
 		  ## mu: features x genotype x allele
346 346
 		  ## Sigma: features x genotype x covariance
347 347
 		  stopifnot(all(copyNumber %in% 0:4))
348
+		  getNu <- function(object){
349
+			  nu[, 1, ] <- nuA(object)
350
+			  nu[, 2, ] <- nuB(object)
351
+			  nu
352
+		  }
353
+		  getPhi <- function(object){
354
+			  phi[,1,] <- phiA(object)
355
+			  phi[,2,] <- phiB(object)
356
+			  phi
357
+		  }
348 358
 		  gts <- lapply(as.list(copyNumber), genotypes)
349 359
 		  nms <- unlist(gts)
350 360
 		  res <- vector("list", length(nms))
351 361
 		  ##names(res) <- paste("copyNumber", copyNumber, sep="")
352 362
 		  names(res) <- nms
363
+		  bnms <- batchNames(object)
364
+		  nu <- array(NA, dim=c(nrow(object), 2, length(bnms)))
365
+		  phi <- array(NA, dim=c(nrow(object), 2, length(bnms)))
353 366
 		  nu <- getNu(object)
354 367
 		  phi <- getPhi(object)
355
-		  mus <- matrix(NA, nrow(nu), 2, dimnames=list(NULL, LETTERS[1:2]))
356
-		  Sigma <- matrix(NA, nrow(nu), 3)
368
+		  ##mus <- matrix(NA, nrow(nu), 2, dimnames=list(NULL, LETTERS[1:2]))
369
+		  mus <- array(NA, dim=c(nrow(nu), 2, length(bnms)),
370
+			       dimnames=list(NULL, LETTERS[1:2],
371
+			       bnms))
372
+		  ## Sigma <- matrix(NA, nrow(nu), 3)
373
+		  Sigma <- array(NA, dim=c(nrow(nu), 3, length(bnms)),
374
+				 dimnames=list(NULL, c("varA", "cor", "varB"),
375
+				 bnms))
357 376
 		  bivariateCenter <- function(nu, phi){
358 377
 			  ##  lexical scope for mus, CA, CB
359
-			  mus[,1] <- log2(nu[1] + CA * phi[1])
360
-			  mus[,2] <- log2(nu[2] + CB * phi[2])
378
+			  mus[,1, ] <- log2(nu[, 1, ] + CA * phi[, 1, ])
379
+			  mus[,2, ] <- log2(nu[, 2, ] + CB * phi[, 2, ])
361 380
 			  mus
362 381
 		  }
363 382
 		  for(i in seq_along(copyNumber)){
... ...
@@ -372,12 +391,24 @@ setMethod("predictionRegion", signature(object="CNSet", copyNumber="integer"),
372 391
 				  gt.corr <- genotypeCorrelation(gt)
373 392
 				  nma <- ifelse(CA == 0, "tau2A.BB", "tau2A.AA")
374 393
 				  nmb <- ifelse(CB == 0, "tau2B.AA", "tau2B.BB")
375
-				  Sigma[,1] <- getVar(object, nma)
376
-				  Sigma[,3] <- getVar(object, nmb)
377
-				  Sigma[,2] <- getCor(object, gt.corr)
394
+				  Sigma[, 1, ] <- getVar(object, nma)
395
+				  Sigma[, 3, ] <- getVar(object, nmb)
396
+				  Sigma[, 2, ] <- getCor(object, gt.corr)
378 397
 				  res[[gt]]$mu <- bivariateCenter(nu, phi)
379 398
 				  res[[gt]]$cov <- Sigma
380 399
 			  }
381 400
 		  }
382 401
 		  return(res)
383 402
 	  })
403
+
404
+
405
+setMethod("xyplotcrlmm", signature(x="formula", data="CNSet", predictRegion="list"),
406
+	  function(x, data, predictRegion, ...){
407
+		  df <- data.frame(A=log2(A(data)), B=log2(B(data)), gt=calls(data), gt.conf=confs(data))#, snp=snpId)
408
+		  ##df <- as.data.frame(data)
409
+		  xyplot(x, df, predictRegion=predictRegion, ...)
410
+	  })
411
+setMethod("xyplot", signature(x="formula", data="CNSet"),
412
+	  function(x, data, ...){
413
+		  xyplotcrlmm(x, data, ...)
414
+})