Browse code

The genotype frequency computed for chromosome X previously used only the women. This is fixed to include the men. If no women are present in a batch, only the genotype frequency is computed -- this is true both for SNPs and nonpolymorphic loci on chromosome X

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

Rob Scharp authored on 16/02/2011 15:57:13
Showing1 changed files

... ...
@@ -1398,12 +1398,13 @@ summarizeNps <- function(strata, index.list, object, batchSize,
1398 1398
 	if(is.lds) {physical <- get("physical"); open(object)}
1399 1399
 	if(verbose) message("      Probe stratum ", strata, " of ", length(index.list))
1400 1400
 	index <- index.list[[strata]]
1401
-	if(CHR.X) {
1402
-		sample.index <- which(object$gender==2)
1403
-		batches <- split(sample.index, as.character(batch(object))[sample.index])
1404
-	} else {
1405
-		batches <- split(seq_along(batch(object)), as.character(batch(object)))
1406
-	}
1401
+##	if(CHR.X) {
1402
+##		sample.index <- which(object$gender==2)
1403
+##		batches <- split(sample.index, as.character(batch(object))[sample.index])
1404
+##	} else {
1405
+##		batches <- split(seq_along(batch(object)), as.character(batch(object)))
1406
+##	}
1407
+	batches <- split(seq_along(batch(object)), as.character(batch(object)))
1407 1408
 	batchnames <- batchNames(object)
1408 1409
 	nr <- length(index)
1409 1410
 	nc <- length(batchnames)
... ...
@@ -1418,15 +1419,16 @@ summarizeNps <- function(strata, index.list, object, batchSize,
1418 1419
 		rm(AVG, BB)
1419 1420
 	}
1420 1421
 	for(k in seq_along(batches)){
1421
-		B <- batches[[k]]
1422
-		N.AA[, k] <- length(B)
1423
-		this.batch <- unique(as.character(batch(object)[B]))
1422
+		sample.index <- batches[[k]]
1423
+		N.AA[, k] <- length(sample.index)
1424
+		if(CHR.X){
1425
+			gender <- object$gender[sample.index]
1426
+			sample.index <- sample.index[gender == 2]
1427
+			if(length(sample.index) == 0) next()
1428
+s		}
1429
+		this.batch <- unique(as.character(batch(object)[sample.index]))
1424 1430
 		j <- match(this.batch, batchnames)
1425
-		I.A <- AA[, B]
1426
-##		if(is.illumina){
1427
-##			I.B <- BB[, B]
1428
-##			I.A <- I.A + I.B
1429
-##		}
1431
+		I.A <- AA[, sample.index]
1430 1432
 		medianA.AA[, k] <- rowMedians(I.A, na.rm=TRUE)
1431 1433
 		madA.AA[, k] <- rowMAD(I.A, na.rm=TRUE)
1432 1434
 		## log2 Transform Intensities
... ...
@@ -1457,9 +1459,10 @@ summarizeSnps <- function(strata,
1457 1459
 ##	} else {
1458 1460
 ##		batches <- split(seq_along(batch(object)), as.character(batch(object)))
1459 1461
 ##	}
1460
-	if(CHR.X){
1461
-		if(verbose) message("        biallelic cluster medians are estimated using only the women for SNPs on chr. X")
1462
-	}
1462
+	## this message can be confusing if no women are in the dataset
1463
+##	if(CHR.X){
1464
+##		if(verbose) message("        biallelic cluster medians are estimated using only the women for SNPs on chr. X")
1465
+##	}
1463 1466
 	batches <- split(seq_along(batch(object)), as.character(batch(object)))
1464 1467
 	batchnames <- batchNames(object)
1465 1468
 	nr <- length(index)