Browse code

bug fig in calculateRTheta

theta.aa, theta.ab, theta.bb should have columns equal to the number of samples in the batch

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

Rob Scharp authored on 01/11/2011 16:31:57
Showing 3 changed files

... ...
@@ -63,7 +63,8 @@ getFeatureData <- function(cdfName, copynumber=FALSE){
63 63
 		M[index, "isSnp"] <- 0L
64 64
 	}
65 65
 	##A few of the snpProbes do not match -- I think it is chromosome Y.
66
-	M[is.na(M[, "chromosome"]), "isSnp"] <- 1L
66
+	if(any(is.na(M[, "chromosome"])))
67
+		M[is.na(M[, "chromosome"]), "isSnp"] <- 1L
67 68
 	M <- data.frame(M)
68 69
 	return(new("AnnotatedDataFrame", data=M))
69 70
 }
... ...
@@ -1383,9 +1384,6 @@ imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
1383 1384
 }
1384 1385
 
1385 1386
 
1386
-
1387
-
1388
-
1389 1387
 ## constructors for Illumina platform
1390 1388
 constructIlluminaFeatureData <- function(gns, cdfName){
1391 1389
 	pkgname <- paste(cdfName, "Crlmm", sep="")
... ...
@@ -1406,6 +1404,8 @@ constructIlluminaFeatureData <- function(gns, cdfName){
1406 1404
 	    data=data.frame(mapping),
1407 1405
 	    varMetadata=data.frame(labelDescription=colnames(mapping)))
1408 1406
 }
1407
+
1408
+
1409 1409
 constructIlluminaAssayData <- function(np, snp, object, storage.mode="environment", nr){
1410 1410
 	stopifnot(identical(snp$gns, featureNames(object)))
1411 1411
 	gns <- c(featureNames(object), np$gns)
... ...
@@ -1669,7 +1669,8 @@ summarizeSnps <- function(strata,
1669 1669
 			  index.list,
1670 1670
 			  object,
1671 1671
 			  GT.CONF.THR,
1672
-			  verbose, CHR.X, ...){
1672
+			  verbose,
1673
+			  CHR.X, ...){
1673 1674
 ##	if(is.lds) {
1674 1675
 ##		physical <- get("physical")
1675 1676
 ##		open(object)
... ...
@@ -2776,6 +2777,9 @@ calculateRTheta <- function(cnSet, genotype=c("AA", "AB", "BB"), batch.name){
2776 2777
 	x[x < 64] <- 64
2777 2778
 	y[y < 64] <- 64
2778 2779
 	theta[, "theta"] <- atan2(y, x)*2/pi
2780
+	if(any(is.na(theta))){
2781
+		stop("NA's in thetas.  Can not compute theta / R values for batches with fewer than 10 samples")
2782
+	}
2779 2783
 	## so that 'R' is available for NP probes
2780 2784
 	y[is.na(y)] <- 0
2781 2785
 	theta[, "R"] <- x+y
... ...
@@ -540,11 +540,12 @@ setMethod("calculateRBaf", signature(object="CNSet"),
540 540
 	RTheta.ab <- calculateRTheta(object, "AB", batch.name)
541 541
 	RTheta.bb <- calculateRTheta(object, "BB", batch.name)
542 542
 
543
-	theta.aa <- matrix(RTheta.aa[, "theta"], nrow(object), ncol(object), byrow=FALSE)
544
-	theta.ab <- matrix(RTheta.ab[, "theta"], nrow(object), ncol(object), byrow=FALSE)
545
-	theta.bb <- matrix(RTheta.bb[, "theta"], nrow(object), ncol(object), byrow=FALSE)
546
-
547 543
 	J <- which(batch(object) == batch.name)
544
+
545
+	theta.aa <- matrix(RTheta.aa[, "theta"], nrow(object), length(J), byrow=FALSE)
546
+	theta.ab <- matrix(RTheta.ab[, "theta"], nrow(object), length(J), byrow=FALSE)
547
+	theta.bb <- matrix(RTheta.bb[, "theta"], nrow(object), length(J), byrow=FALSE)
548
+
548 549
 	a <- A(object)[, J, drop=FALSE]
549 550
 	b <- B(object)[, J, drop=FALSE] ## NA's for b where nonpolymorphic
550 551
 	dns <- dimnames(a)
... ...
@@ -74,7 +74,8 @@ pathToCels <- "/thumper/ctsa/snpmicroarray/hapmap/raw/affy/1m"
74 74
 if(getRversion() < "2.13.0"){
75 75
 	rpath <- getRversion()
76 76
 } else rpath <- "trunk"
77
-outdir <- paste("/amber1/archive/oncbb/rscharpf/crlmm_vignette/", rpath, sep="")
77
+##outdir <- paste("/amber1/archive/oncbb/rscharpf/crlmm_vignette/", rpath, sep="")
78
+outdir <- paste("/amber1/archive/oncbb/rscharpf/crlmm_vignette/", rpath, "/gw6v1.5/", sep="")
78 79
 dir.create(outdir, recursive=TRUE, showWarnings=FALSE)
79 80
 @
80 81