Browse code

update calculateRTheta for NP markers

use atan2(x,x)*2/pi=1

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

Rob Scharp authored on 24/11/2011 12:50:18
Showing2 changed files

... ...
@@ -2744,21 +2744,23 @@ ABpanel <- function(x, y, predictRegion,
2744 2744
 		panel.xyplot(x, y, ...)
2745 2745
 	}
2746 2746
 	pn <- panel.number()
2747
-	for(CN in copyNumber){
2748
-		gts <- genotypes(CN)
2749
-		index <- match(gts, names(predictRegion))
2750
-		pr <- predictRegion[index]
2751
-		for(i in seq_along(pr)){
2752
-			## scale?
2753
-			pr2 <- pr[[i]]
2754
-			mu <- pr2$mu[pn, , , drop=FALSE] ## pn=panel number
2755
-			Sigma <- pr2$cov[pn, , ,drop=FALSE]
2756
-			for(j in seq_len(dim(mu)[3])){
2757
-				dat.ellipse <- ellipse(x=Sigma[, 2, j],
2758
-						       centre=mu[ , , j],
2759
-						       scale=c(sqrt(Sigma[ , 1, j]),
2760
-						       sqrt(Sigma[, 3, j])))
2761
-				lpolygon(dat.ellipse[,1], dat.ellipse[,2], ...)
2747
+	if(!is.null(copyNumber)){
2748
+		for(CN in copyNumber){
2749
+			gts <- genotypes(CN)
2750
+			index <- match(gts, names(predictRegion))
2751
+			pr <- predictRegion[index]
2752
+			for(i in seq_along(pr)){ ## for each genotype
2753
+				## scale?
2754
+				pr2 <- pr[[i]]
2755
+				mu <- pr2$mu[pn, , , drop=FALSE] ## pn=panel number
2756
+				Sigma <- pr2$cov[pn, , ,drop=FALSE]
2757
+				for(j in seq_len(dim(mu)[3])){  ## for each batch
2758
+					dat.ellipse <- ellipse(x=Sigma[, 2, j],
2759
+							       centre=mu[ , , j],
2760
+							       scale=c(sqrt(Sigma[ , 1, j]),
2761
+							       sqrt(Sigma[, 3, j])))
2762
+					lpolygon(dat.ellipse[,1], dat.ellipse[,2], ...)
2763
+				}
2762 2764
 			}
2763 2765
 		}
2764 2766
 	}
... ...
@@ -2777,12 +2779,16 @@ calculateRTheta <- function(cnSet, genotype=c("AA", "AB", "BB"), batch.name){
2777 2779
 	x[x < 64] <- 64
2778 2780
 	y[y < 64] <- 64
2779 2781
 	theta[, "theta"] <- atan2(y, x)*2/pi
2780
-	if(any(is.na(theta[, "theta"]))){
2781
-		stop("NA's in thetas.  Can not compute theta / R values for batches with fewer than 10 samples")
2782
-	}
2782
+##	if(any(is.na(theta[, "theta"]))){
2783
+##		##stop("NA's in thetas.  Can not compute theta / R values for batches with fewer than 10 samples")
2784
+##	}
2783 2785
 	## so that 'R' is available for NP probes
2784 2786
 	y[is.na(y)] <- 0
2785 2787
 	theta[, "R"] <- x+y
2788
+	theta[!isSnp(cnSet), 1] <- 1
2789
+	if(any(is.na(theta[, "theta"]))){
2790
+		stop("NA's in thetas.  Can not compute theta / R values for batches with fewer than 10 samples")
2791
+	}
2786 2792
 	return(theta)
2787 2793
 }
2788 2794
 
... ...
@@ -519,6 +519,7 @@ setMethod("xyplotcrlmm", signature(x="formula", data="CNSet", predictRegion="lis
519 519
 		  ##predictRegion is an argument of ABpanel
520 520
 		  xyplot(x, df, predictRegion=predictRegion, ...)
521 521
 	  })
522
+
522 523
 setMethod("xyplot", signature(x="formula", data="CNSet"),
523 524
 	  function(x, data, ...){
524 525
 		  if("predictRegion" %in% names(list(...))){
... ...
@@ -552,6 +553,8 @@ setMethod("calculateRBaf", signature(object="CNSet"),
552 553
 
553 554
 		  a <- A(object)[, J, drop=FALSE]
554 555
 		  b <- B(object)[, J, drop=FALSE] ## NA's for b where nonpolymorphic
556
+		  is.np <- !isSnp(object)
557
+		  b[is.np, ] <- a[is.np, ]
555 558
 		  dns <- dimnames(a)
556 559
 		  dimnames(a) <- dimnames(b) <- NULL
557 560
 		  obs.theta <- atan2(b, a)*2/pi
... ...
@@ -588,7 +591,7 @@ setMethod("calculateRBaf", signature(object="CNSet"),
588 591
 		  r.expected[lessAA] <- r.aa[lessAA]
589 592
 		  r.expected[grBB] <- r.bb[grBB]
590 593
 
591
-		  index.np <- which(!isSnp(object))
594
+		  index.np <- which(is.np)
592 595
 		  a.np <- A(object)[index.np, ]
593 596
 		  meds <- rowMedians(a.np, na.rm=TRUE)
594 597
 		  r.expected[index.np, ] <- matrix(meds, length(index.np), ncol(a.np))