Browse code

Add calculateRBaf method

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

Rob Scharp authored on 01/10/2011 04:50:37
Showing5 changed files

... ...
@@ -1,7 +1,7 @@
1 1
 Package: crlmm
2 2
 Type: Package
3 3
 Title: Genotype Calling (CRLMM) and Copy Number Analysis tool for Affymetrix SNP 5.0 and 6.0 and Illumina arrays.
4
-Version: 1.11.41
4
+Version: 1.11.42
5 5
 Date: 2010-12-10
6 6
 Author: Benilton S Carvalho <Benilton.Carvalho@cancer.org.uk>, Robert Scharpf <rscharpf@jhsph.edu>, Matt Ritchie <mritchie@wehi.edu.au>, Ingo Ruczinski <iruczins@jhsph.edu>, Rafael A Irizarry
7 7
 Maintainer: Benilton S Carvalho <Benilton.Carvalho@cancer.org.uk>, Robert Scharpf <rscharpf@jhsph.edu>, Matt Ritchie <mritchie@wehi.EDU.AU>
... ...
@@ -64,6 +64,7 @@ export(crlmm,
64 64
        genotype.Illumina,
65 65
        crlmmCopynumber2, crlmmCopynumberLD, crlmmCopynumber)
66 66
 export(genotypes, totalCopynumber, rawCopynumber, xyplot)
67
-exportMethods(A, B, calculatePosteriorMean, corr, nuA, nuB, phiA, phiB, predictionRegion, posteriorProbability, tau2, Ns, medians, mads, xyplot)
67
+exportMethods(A, B, calculatePosteriorMean, corr, nuA, nuB, phiA, phiB, predictionRegion, posteriorProbability, tau2, Ns, medians, mads,
68
+	      xyplot, calculateRBaf)
68 69
 export(ABpanel, constructInf, preprocessInf, genotypeInf)
69 70
 exportClasses(PredictionRegion)
... ...
@@ -104,3 +104,4 @@ setGeneric("predictionRegion", function(object, copyNumber=0:4)
104 104
 	   standardGeneric("predictionRegion"))
105 105
 setGeneric("xyplot", useAsDefault=function(x, data, ...) lattice::xyplot(x, data,...))
106 106
 setGeneric("xyplotcrlmm", function(x, data, predictRegion,...) standardGeneric("xyplotcrlmm"))
107
+setGeneric("calculateRBaf", function(object, batch.name) standardGeneric("calculateRBaf"))
... ...
@@ -1306,6 +1306,8 @@ imputeCenter <- function(muA, muB, index.complete, index.missing){
1306 1306
 			betahat <- solve(crossprod(X), crossprod(X,Y))
1307 1307
 			X <- cbind(1, mnA[index[[j]],  -j, drop=FALSE],  mnB[index[[j]],  -j, drop=FALSE])
1308 1308
 			mus <- X %*% betahat
1309
+			## threshold imputed intensities that are very small
1310
+			mus[mus < 64] <- 64
1309 1311
 			muA[index[[j]], j] <- mus[, 1]
1310 1312
 			muB[index[[j]], j] <- mus[, 2]
1311 1313
 		}
... ...
@@ -2764,6 +2766,25 @@ ABpanel <- function(x, y, predictRegion,
2764 2766
 	}
2765 2767
 }
2766 2768
 
2769
+calculateRTheta <- function(cnSet, genotype=c("AA", "AB", "BB"), batch.name){
2770
+	genotype <- match.arg(genotype)
2771
+	j <- match(batch.name, batchNames(cnSet))
2772
+	centers <- medians(cnSet, i=snp.index, j)
2773
+	theta <- matrix(NA, nrow(cnSet), 2)
2774
+	colnames(theta) <- c("theta", "R")
2775
+	x <- centers[, "A", genotype, 1]
2776
+	y <- centers[, "B", genotype, 1]
2777
+	## small imputed values -- should fix imputeCenter
2778
+	x[x < 64] <- 64
2779
+	y[y < 64] <- 64
2780
+	theta[, "theta"] <- atan2(y, x)*2/pi
2781
+	theta[, "R"] <- x+y
2782
+	return(theta)
2783
+}
2784
+
2785
+
2786
+
2787
+
2767 2788
 calculateBAFandLRR <- function(cnSet) {
2768 2789
   ## Calculates B allele frequency and log2 R ratio for all snps for a given cnSet
2769 2790
   ## Returns a 3-D array of size Num Of Snps x Num of Samples x 2
... ...
@@ -528,3 +528,69 @@ setMethod("xyplot", signature(x="formula", data="CNSet"),
528 528
 			  callNextMethod()
529 529
 		  }
530 530
 })
531
+
532
+setMethod("calculateRBaf", signature(object="CNSet"),
533
+	  function(object, batch.name){
534
+	if(missing(batch.name)) batch.name <- batchNames(object)[1]
535
+	stopifnot(batch.name %in% batchNames(object))
536
+	if(length(batch.name) > 1){
537
+		warning("only the first batch in batch.name processed")
538
+		batch.name <- batch.name[1]
539
+	}
540
+	RTheta.aa <- calculateRTheta(object, "AA", batch.name)
541
+	RTheta.ab <- calculateRTheta(object, "AB", batch.name)
542
+	RTheta.bb <- calculateRTheta(object, "BB", batch.name)
543
+
544
+	theta.aa <- matrix(RTheta.aa[, "theta"], nrow(object), ncol(object), byrow=FALSE)
545
+	theta.ab <- matrix(RTheta.ab[, "theta"], nrow(object), ncol(object), byrow=FALSE)
546
+	theta.bb <- matrix(RTheta.bb[, "theta"], nrow(object), ncol(object), byrow=FALSE)
547
+
548
+	J <- which(batch(object) == batch.name)
549
+	a <- A(object)[, J, drop=FALSE]
550
+	b <- B(object)[, J, drop=FALSE] ## NA's for b where nonpolymorphic
551
+	dimnames(a) <- dimnames(b) <- NULL
552
+	obs.theta <- atan2(b, a)*2/pi
553
+
554
+	lessAA <- obs.theta < theta.aa
555
+	lessAB <- obs.theta < theta.ab
556
+	lessBB <- obs.theta < theta.bb
557
+	grAA <- !lessAA
558
+	grAB <- !lessAB
559
+	grBB <- !lessBB
560
+	not.na <- !is.na(theta.aa)
561
+	I1 <- grAA & lessAB & not.na
562
+	I2 <- grAB & lessBB & not.na
563
+
564
+	bf <- matrix(NA, nrow(object), ncol(a))
565
+	bf[I1] <- 0.5 * ((obs.theta-theta.aa)/(theta.ab-theta.aa))[I1]
566
+	bf[I2] <- (.5 * (obs.theta - theta.ab) / (theta.bb - theta.ab))[I2] + 0.5
567
+	bf[lessAA] <- 0
568
+	bf[grBB] <- 1
569
+
570
+	r.expected <- matrix(NA, nrow(object), ncol(a))
571
+	r.aa <- matrix(RTheta.aa[, "R"], nrow(object), ncol(object), byrow=FALSE)
572
+	r.ab <- matrix(RTheta.ab[, "R"], nrow(object), ncol(object), byrow=FALSE)
573
+	r.bb <- matrix(RTheta.bb[, "R"], nrow(object), ncol(object), byrow=FALSE)
574
+	rm(RTheta.aa, RTheta.ab, RTheta.bb); gc()
575
+	obs.r <- a+b
576
+
577
+	lessAA <- lessAA & not.na
578
+	grBB <- grBB & not.na
579
+	tmp <- ((obs.theta - theta.aa) * (r.ab-r.aa)/(theta.ab-theta.aa))[I1] + r.aa[I1]
580
+	r.expected[I1] <- tmp
581
+	tmp <- ((obs.theta - theta.ab) * (r.bb - r.ab)/(theta.bb-theta.ab))[I2] + r.ab[I2]
582
+	r.expected[I2] <- tmp
583
+	r.expected[lessAA] <- r.aa[lessAA]
584
+	r.expected[grBB] <- r.bb[grBB]
585
+
586
+	index.np <- which(!isSnp(object))
587
+	a.np <- A(object)[index.np, ]
588
+	meds <- rowMedians(a.np, na.rm=TRUE)
589
+	r.expected[index.np, ] <- matrix(meds, length(index.np), ncol(a.np))
590
+	lrr <- log2(obs.r/r.expected)
591
+
592
+	dimnames(bf) <- dimnames(lrr) <- dimnames(a)
593
+	res <- list(baf=bf,
594
+		    lrr=lrr)
595
+	return(res)
596
+})