Browse code

Check that only autosomes are in object passed to calculateRBaf

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

Rob Scharp authored on 01/11/2011 16:32:14
Showing2 changed files

... ...
@@ -2771,8 +2771,8 @@ calculateRTheta <- function(cnSet, genotype=c("AA", "AB", "BB"), batch.name){
2771 2771
 	centers <- medians(cnSet, i=seq_len(nrow(cnSet)), j)
2772 2772
 	theta <- matrix(NA, nrow(cnSet), 2)
2773 2773
 	colnames(theta) <- c("theta", "R")
2774
-	x <- centers[, "A", genotype, 1]
2775
-	y <- centers[, "B", genotype, 1]
2774
+	x <- centers[, "A", genotype, j]
2775
+	y <- centers[, "B", genotype, j]
2776 2776
 	## small imputed values -- should fix imputeCenter
2777 2777
 	x[x < 64] <- 64
2778 2778
 	y[y < 64] <- 64
... ...
@@ -530,68 +530,72 @@ setMethod("xyplot", signature(x="formula", data="CNSet"),
530 530
 
531 531
 setMethod("calculateRBaf", signature(object="CNSet"),
532 532
 	  function(object, batch.name){
533
-	if(missing(batch.name)) batch.name <- batchNames(object)[1]
534
-	stopifnot(batch.name %in% batchNames(object))
535
-	if(length(batch.name) > 1){
536
-		warning("only the first batch in batch.name processed")
537
-		batch.name <- batch.name[1]
538
-	}
539
-	RTheta.aa <- calculateRTheta(object, "AA", batch.name)
540
-	RTheta.ab <- calculateRTheta(object, "AB", batch.name)
541
-	RTheta.bb <- calculateRTheta(object, "BB", batch.name)
533
+		  all.autosomes <- all(chromosome(object) < 23)
534
+		  if(!all.autosomes){
535
+			  stop("method currently only defined for chromosomes 1-22")
536
+		  }
537
+		  if(missing(batch.name)) batch.name <- batchNames(object)[1]
538
+		  stopifnot(batch.name %in% batchNames(object))
539
+		  if(length(batch.name) > 1){
540
+			  warning("only the first batch in batch.name processed")
541
+			  batch.name <- batch.name[1]
542
+		  }
543
+		  RTheta.aa <- calculateRTheta(object, "AA", batch.name)
544
+		  RTheta.ab <- calculateRTheta(object, "AB", batch.name)
545
+		  RTheta.bb <- calculateRTheta(object, "BB", batch.name)
542 546
 
543
-	J <- which(batch(object) == batch.name)
547
+		  J <- which(batch(object) == batch.name)
544 548
 
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)
549
+		  theta.aa <- matrix(RTheta.aa[, "theta"], nrow(object), length(J), byrow=FALSE)
550
+		  theta.ab <- matrix(RTheta.ab[, "theta"], nrow(object), length(J), byrow=FALSE)
551
+		  theta.bb <- matrix(RTheta.bb[, "theta"], nrow(object), length(J), byrow=FALSE)
548 552
 
549
-	a <- A(object)[, J, drop=FALSE]
550
-	b <- B(object)[, J, drop=FALSE] ## NA's for b where nonpolymorphic
551
-	dns <- dimnames(a)
552
-	dimnames(a) <- dimnames(b) <- NULL
553
-	obs.theta <- atan2(b, a)*2/pi
553
+		  a <- A(object)[, J, drop=FALSE]
554
+		  b <- B(object)[, J, drop=FALSE] ## NA's for b where nonpolymorphic
555
+		  dns <- dimnames(a)
556
+		  dimnames(a) <- dimnames(b) <- NULL
557
+		  obs.theta <- atan2(b, a)*2/pi
554 558
 
555
-	lessAA <- obs.theta < theta.aa
556
-	lessAB <- obs.theta < theta.ab
557
-	lessBB <- obs.theta < theta.bb
558
-	grAA <- !lessAA
559
-	grAB <- !lessAB
560
-	grBB <- !lessBB
561
-	not.na <- !is.na(theta.aa)
562
-	I1 <- grAA & lessAB & not.na
563
-	I2 <- grAB & lessBB & not.na
559
+		  lessAA <- obs.theta < theta.aa
560
+		  lessAB <- obs.theta < theta.ab
561
+		  lessBB <- obs.theta < theta.bb
562
+		  grAA <- !lessAA
563
+		  grAB <- !lessAB
564
+		  grBB <- !lessBB
565
+		  not.na <- !is.na(theta.aa)
566
+		  I1 <- grAA & lessAB & not.na
567
+		  I2 <- grAB & lessBB & not.na
564 568
 
565
-	bf <- matrix(NA, nrow(object), ncol(a))
566
-	bf[I1] <- 0.5 * ((obs.theta-theta.aa)/(theta.ab-theta.aa))[I1]
567
-	bf[I2] <- (.5 * (obs.theta - theta.ab) / (theta.bb - theta.ab))[I2] + 0.5
568
-	bf[lessAA] <- 0
569
-	bf[grBB] <- 1
569
+		  bf <- matrix(NA, nrow(object), ncol(a))
570
+		  bf[I1] <- 0.5 * ((obs.theta-theta.aa)/(theta.ab-theta.aa))[I1]
571
+		  bf[I2] <- (.5 * (obs.theta - theta.ab) / (theta.bb - theta.ab))[I2] + 0.5
572
+		  bf[lessAA] <- 0
573
+		  bf[grBB] <- 1
570 574
 
571
-	r.expected <- matrix(NA, nrow(object), ncol(a))
572
-	r.aa <- matrix(RTheta.aa[, "R"], nrow(object), ncol(object), byrow=FALSE)
573
-	r.ab <- matrix(RTheta.ab[, "R"], nrow(object), ncol(object), byrow=FALSE)
574
-	r.bb <- matrix(RTheta.bb[, "R"], nrow(object), ncol(object), byrow=FALSE)
575
-	rm(RTheta.aa, RTheta.ab, RTheta.bb); gc()
576
-	obs.r <- a+b
575
+		  r.expected <- matrix(NA, nrow(object), ncol(a))
576
+		  r.aa <- matrix(RTheta.aa[, "R"], nrow(object), ncol(object), byrow=FALSE)
577
+		  r.ab <- matrix(RTheta.ab[, "R"], nrow(object), ncol(object), byrow=FALSE)
578
+		  r.bb <- matrix(RTheta.bb[, "R"], nrow(object), ncol(object), byrow=FALSE)
579
+		  rm(RTheta.aa, RTheta.ab, RTheta.bb); gc()
580
+		  obs.r <- a+b
577 581
 
578
-	lessAA <- lessAA & not.na
579
-	grBB <- grBB & not.na
580
-	tmp <- ((obs.theta - theta.aa) * (r.ab-r.aa)/(theta.ab-theta.aa))[I1] + r.aa[I1]
581
-	r.expected[I1] <- tmp
582
-	tmp <- ((obs.theta - theta.ab) * (r.bb - r.ab)/(theta.bb-theta.ab))[I2] + r.ab[I2]
583
-	r.expected[I2] <- tmp
584
-	r.expected[lessAA] <- r.aa[lessAA]
585
-	r.expected[grBB] <- r.bb[grBB]
582
+		  lessAA <- lessAA & not.na
583
+		  grBB <- grBB & not.na
584
+		  tmp <- ((obs.theta - theta.aa) * (r.ab-r.aa)/(theta.ab-theta.aa))[I1] + r.aa[I1]
585
+		  r.expected[I1] <- tmp
586
+		  tmp <- ((obs.theta - theta.ab) * (r.bb - r.ab)/(theta.bb-theta.ab))[I2] + r.ab[I2]
587
+		  r.expected[I2] <- tmp
588
+		  r.expected[lessAA] <- r.aa[lessAA]
589
+		  r.expected[grBB] <- r.bb[grBB]
586 590
 
587
-	index.np <- which(!isSnp(object))
588
-	a.np <- A(object)[index.np, ]
589
-	meds <- rowMedians(a.np, na.rm=TRUE)
590
-	r.expected[index.np, ] <- matrix(meds, length(index.np), ncol(a.np))
591
-	lrr <- log2(obs.r/r.expected)
591
+		  index.np <- which(!isSnp(object))
592
+		  a.np <- A(object)[index.np, ]
593
+		  meds <- rowMedians(a.np, na.rm=TRUE)
594
+		  r.expected[index.np, ] <- matrix(meds, length(index.np), ncol(a.np))
595
+		  lrr <- log2(obs.r/r.expected)
592 596
 
593
-	dimnames(bf) <- dimnames(lrr) <- dns
594
-	res <- list(baf=bf,
595
-		    lrr=lrr)
596
-	return(res)
597
-})
597
+		  dimnames(bf) <- dimnames(lrr) <- dns
598
+		  res <- list(baf=bf,
599
+			      lrr=lrr)
600
+		  return(res)
601
+	  })