Browse code

CN estimation for Chromosome X will not throw an error when too few men or women are available in a batch.

Instead, a message is printed and NAs (likely) will result in the CN estimate

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

Rob Scharp authored on 04/10/2010 02:22:13
Showing1 changed files

... ...
@@ -171,7 +171,7 @@ genotype <- function(filenames,
171 171
 		       seed=seed,
172 172
 		       verbose=verbose)
173 173
 	if(verbose) message("Saving callSet.rda")
174
-	save(callSet, file=file.path(outdir, "callSet.rda"))
174
+	##save(callSet, file=file.path(outdir, "callSet.rda"))
175 175
 	if(!is.lds) A(callSet) <- AA
176 176
 	## otherwise, the normalized values were written to file... no need to do anything
177 177
 	rm(AA)
... ...
@@ -901,12 +901,13 @@ fit.lm4 <- function(strata,
901 901
 		X <- cbind(1, log2(c(tmp[, 1], tmp[, 2])))
902 902
 		Y <- log2(c(tmp[, 3], tmp[, 4]))
903 903
 		betahat <- solve(crossprod(X), crossprod(X, Y))
904
-		X.men <- cbind(1, medianA.AA.M[, k])
905
-		Yhat1 <- as.numeric(X.men %*% betahat)
906
-		## put intercept and slope for men in nuB and phiB
907
-		phiB[, k] <- 2^(Yhat1)
908
-		nuB[, k] <- 2^(medianA.AA.M[, k]) - phiB[, k]
909
-
904
+		if(enough.men){
905
+			X.men <- cbind(1, medianA.AA.M[, k])
906
+			Yhat1 <- as.numeric(X.men %*% betahat)
907
+			## put intercept and slope for men in nuB and phiB
908
+			phiB[, k] <- 2^(Yhat1)
909
+			nuB[, k] <- 2^(medianA.AA.M[, k]) - phiB[, k]
910
+		}
910 911
 		X.fem <- cbind(1, medianA.AA.F[, k])
911 912
 		Yhat2 <- as.numeric(X.fem %*% betahat)
912 913
 		phiA[, k] <- 2^(Yhat2)
... ...
@@ -1492,7 +1493,8 @@ shrinkSummary <- function(object,
1492 1493
 	if(type[[1]] == "X.SNP"){
1493 1494
 		gender <- object$gender
1494 1495
 		if(sum(gender == 2) < 3) {
1495
-			return("too few females to estimate within genotype summary statistics on CHR X")
1496
+			message("too few females to estimate within genotype summary statistics on CHR X")
1497
+			return(object)
1496 1498
 		}
1497 1499
 		CHR.X <- TRUE
1498 1500
 	} else CHR.X <- FALSE
... ...
@@ -1540,7 +1542,8 @@ genotypeSummary <- function(object,
1540 1542
 	if(type == "X.SNP" | type=="X.NP"){
1541 1543
 		gender <- object$gender
1542 1544
 		if(sum(gender == 2) < 3) {
1543
-			return("too few females to estimate within genotype summary statistics on CHR X")
1545
+			message("too few females to estimate within genotype summary statistics on CHR X")
1546
+			return(object)
1544 1547
 		}
1545 1548
 		CHR.X <- TRUE
1546 1549
 	} else CHR.X <- FALSE