Browse code

Update posteriorProbability for handling nonpolymorphic markers

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

Rob Scharp authored on 01/10/2011 04:49:02
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.28
4
+Version: 1.11.29
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>
... ...
@@ -298,36 +298,69 @@ setMethod("posteriorProbability", signature(object="CNSet"),
298 298
 		  prob <- array(NA, dim=c(nrow(object), ncol(object), length(copyNumber)),
299 299
 				dimnames=list(NULL,NULL, paste("copynumber",copyNumber, sep="")))
300 300
 		  bns <- batchNames(object)
301
-		  for(i in seq_along(copyNumber)){
302
-			  G <- genotypes(copyNumber[i])
303
-			  P <- array(NA, dim=c(nrow(object), ncol(object), length(G)),
304
-				     dimnames=list(NULL,NULL,G))
305
-			  for(g in seq_along(G)){
306
-				  gt <- G[g]
307
-				  for(j in seq_along(bns)){
308
-					  this.batch <- bns[j]
309
-					  sample.index <- which(batch(object) == this.batch)
310
-					  mu <- predictRegion[[gt]]$mu[, , this.batch, drop=FALSE]
311
-					  dim(mu) <- dim(mu)[1:2]
312
-					  if(all(is.na(mu))){
313
-						  P[, sample.index, gt] <- NA
314
-					  } else {
315
-						  Sigma <- predictRegion[[gt]]$cov[, , this.batch, drop=FALSE]
316
-						  dim(Sigma) <- dim(Sigma)[1:2]
317
-						  P[, sample.index, gt] <- dbvn(x=logI[, sample.index, , drop=FALSE], mu=mu, Sigma=Sigma)
301
+		  snp.index <- which(isSnp(object))
302
+		  if(length(snp.index) > 0){
303
+			  for(i in seq_along(copyNumber)){
304
+				  G <- genotypes(copyNumber[i])
305
+				  P <- array(NA, dim=c(length(snp.index), ncol(object), length(G)),
306
+					     dimnames=list(NULL,NULL,G))
307
+				  for(g in seq_along(G)){
308
+					  gt <- G[g]
309
+					  for(j in seq_along(bns)){
310
+						  this.batch <- bns[j]
311
+						  sample.index <- which(batch(object) == this.batch)
312
+						  mu <- predictRegion[[gt]]$mu[snp.index, , this.batch, drop=FALSE]
313
+						  dim(mu) <- dim(mu)[1:2]
314
+						  if(all(is.na(mu))){
315
+							  P[, sample.index, gt] <- NA
316
+						  } else {
317
+							  Sigma <- predictRegion[[gt]]$cov[snp.index, , this.batch, drop=FALSE]
318
+							  dim(Sigma) <- dim(Sigma)[1:2]
319
+							  P[, sample.index, gt] <- dbvn(x=logI[snp.index, sample.index, , drop=FALSE], mu=mu, Sigma=Sigma)
320
+						  }
318 321
 					  }
319 322
 				  }
323
+				  ## the marginal probability for the total copy number
324
+				  ##  -- integrate out the probability for the different genotypes
325
+				  for(j in 1:ncol(P)){
326
+					  PP <- P[, j, , drop=FALSE]
327
+					  dim(PP) <- dim(PP)[c(1, 3)]
328
+					  prob[snp.index, j, i] <- rowSums(PP, na.rm=TRUE)
329
+				  }
320 330
 			  }
321
-			  ## the marginal probability for the total copy number
322
-			  ##  -- integrate out the probability for the different genotypes
323
-			  for(j in 1:ncol(P)){
324
-				  PP <- P[, j, , drop=FALSE]
325
-				  dim(PP) <- dim(PP)[c(1, 3)]
326
-				  prob[, j, i] <- rowSums(PP, na.rm=TRUE)
327
-			  }
328
-		  }
329
-
330
-		  ## scale by prior weights and divide by normalizing
331
+		  } ## length(snp.index) > 0
332
+		  ##---------------------------------------------------------------------------
333
+		  ##
334
+		  ##  nonpolymorphic markers
335
+		  ##
336
+		  ##---------------------------------------------------------------------------
337
+		  np.index <- which(!isSnp(object))
338
+		  if(length(np.index) > 0){
339
+			  for(i in seq_along(copyNumber)){
340
+				  G <- genotypes(copyNumber[i], is.snp=FALSE)
341
+				  P <- array(NA, dim=c(length(np.index), ncol(object), length(G)),
342
+					     dimnames=list(NULL,NULL,G))
343
+				  for(g in seq_along(G)){
344
+					  gt <- G[g]
345
+					  for(j in seq_along(bns)){
346
+						  this.batch <- bns[j]
347
+						  sample.index <- which(batch(object) == this.batch)
348
+						  mu <- predictRegion[[gt]]$mu[np.index, 1, this.batch, drop=FALSE]
349
+						  dim(mu) <- dim(mu)[1]
350
+						  if(all(is.na(mu))){
351
+							  P[, sample.index, gt] <- NA
352
+						  } else {
353
+							  Sigma <- predictRegion[[gt]]$cov[np.index, 1, this.batch, drop=FALSE]
354
+							  dim(Sigma) <- dim(Sigma)[1]
355
+							  P[, sample.index, gt] <- dnorm(logI[np.index, sample.index, 1], mean=mu, sd=Sigma)
356
+						  }
357
+					  } ## j in seq_along(bns)
358
+				  } ## g in seq_along(G)
359
+				  ## the marginal probability for the total copy number
360
+				  ##  -- integrate out the probability for the different genotypes
361
+				  prob[np.index, , i] <- P
362
+			  } ## seq_along(copyNumber)
363
+		  } ## length(np.index) > 0
331 364
 		  ##constant.  Probs across copy number states must
332 365
 		  ##sum to one.  nc <- apply(prob, c(1,3), sum,
333 366
 		  ##na.rm=TRUE)
... ...
@@ -399,6 +432,7 @@ setMethod("predictionRegion", signature(object="CNSet", copyNumber="integer"),
399 432
 			  mus[,2, ] <- log2(nu[, 2, ] + CB * phi[, 2, ])
400 433
 			  mus
401 434
 		  }
435
+		  np.index <- which(!isSnp(object))
402 436
 		  for(i in seq_along(copyNumber)){
403 437
 			  G <- genotypes(copyNumber[i])
404 438
 			  tmp <- vector("list", length(G))
... ...
@@ -414,6 +448,10 @@ setMethod("predictionRegion", signature(object="CNSet", copyNumber="integer"),
414 448
 				  Sigma[, 1, ] <- getVar(object, nma)
415 449
 				  Sigma[, 3, ] <- getVar(object, nmb)
416 450
 				  Sigma[, 2, ] <- getCor(object, gt.corr)
451
+				  if(length(np.index) > 0){
452
+					  Sigma[, 1, ] <- getVar(object, "tau2A.AA")
453
+					  Sigma[np.index, 2, ] <- NA
454
+				  }
417 455
 				  res[[gt]]$mu <- bivariateCenter(nu, phi)
418 456
 				  res[[gt]]$cov <- Sigma
419 457
 			  }