Browse code

Add generics and methods for predictionRegion, posteriorProbability, and calculatePosteriorMean.

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

Rob Scharp authored on 01/10/2011 04:46:51
Showing6 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.9
4
+Version: 1.11.10
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>
... ...
@@ -34,6 +34,7 @@ Collate: AllGenerics.R
34 34
 	 methods-eSet.R
35 35
          methods-SnpSuperSet.R
36 36
          cnrma-functions.R
37
+	 cnset-accessors.R
37 38
          crlmm-functions.R
38 39
 	 crlmmGT2.R
39 40
          crlmm-illumina.R
... ...
@@ -74,6 +74,6 @@ export(crlmm,
74 74
        genotype2, genotypeLD,
75 75
        genotype.Illumina,
76 76
        crlmmCopynumber2, crlmmCopynumberLD, crlmmCopynumber)
77
-export(totalCopynumber, rawCopynumber)
78
-exportMethods(A, B, nuA, nuB, phiA, phiB, corr, tau2, Ns, medians, mads)
77
+export(genotypes, totalCopynumber, rawCopynumber)
78
+exportMethods(A, B, calculatePosteriorMean, corr, nuA, nuB, phiA, phiB, predictionRegion, posteriorProbability, tau2, Ns, medians, mads)
79 79
 export(constructInf, preprocessInf, genotypeInf)
... ...
@@ -95,5 +95,10 @@ setGeneric("flags<-", function(object, value) standardGeneric("flags<-"))
95 95
 setGeneric("posteriorMean", function(object) standardGeneric("posteriorMean"))
96 96
 setGeneric("posteriorMean<-", function(object, value) standardGeneric("posteriorMean<-"))
97 97
 
98
+setGeneric("posteriorProbability", function(object, predictRegion, copyNumber=0:4)
99
+	   standardGeneric("posteriorProbability"))
100
+setGeneric("calculatePosteriorMean", function(object, posteriorProb, copyNumber=0:4, w, ...)
101
+	   standardGeneric("calculatePosteriorMean"))
98 102
 
99
-
103
+setGeneric("predictionRegion", function(object, copyNumber=0:4)
104
+	   standardGeneric("predictionRegion"))
... ...
@@ -2346,79 +2346,79 @@ imputeAcrossBatch <- function(N.AA, N.AB, N.BB,
2346 2346
 }
2347 2347
 
2348 2348
 
2349
-calculatePosteriorMean <- function(object, type=c("SNP", "NP", "X.SNP", "X.NP"), verbose=TRUE,
2350
-				   prior.prob=c(1/7, 1/7, 3/7, 1/7, 1/7),
2351
-				   CN=0:4, scale.sd=1){
2352
-	stopifnot(type %in% c("SNP", "NP", "X.SNP", "X.NP"))
2353
-	stopifnot(sum(prior.prob)==1)
2354
-	stopifnot(length(CN) == length(prior.prob))
2355
-	batch <- batch(object)
2356
-	is.snp <- isSnp(object)
2357
-	is.autosome <- chromosome(object) < 23
2358
-	is.annotated <- !is.na(chromosome(object))
2359
-	is.X <- chromosome(object) == 23
2360
-	is.lds <- is(calls(object), "ffdf") | is(calls(object), "ff_matrix")
2361
-	if(is.lds) stopifnot(isPackageLoaded("ff"))
2362
-	samplesPerBatch <- table(as.character(batch(object)))
2363
-	if(!"posteriorMean" %in% assayDataElementNames(object)){
2364
-		message("adding <posteriorMean> slot to assayData")
2365
-		pM <- matrix(NA, nrow(object), ncol(object), dimnames=list(featureNames(object), sampleNames(object)))
2366
-		tmp <- assayDataNew(alleleA=A(object),
2367
-				    alleleB=B(object),
2368
-				    call=calls(object),
2369
-				    callProbability=snpCallProbability(object),
2370
-				    posteriorMean=pM)
2371
-		assayData(object) <- tmp
2372
-	}
2373
-	## add new assay data element for posterior probabilities
2374
-	mylabel <- function(marker.type){
2375
-		switch(marker.type,
2376
-		       SNP="autosomal SNPs",
2377
-		       NP="autosomal nonpolymorphic markers",
2378
-		       X.SNP="chromosome X SNPs",
2379
-		       X.NP="chromosome X nonpolymorphic markers")
2380
-	}
2381
-	if(type=="SNP"){
2382
-		if(verbose) message(paste("...", mylabel(type)))
2383
-		marker.index <- whichMarkers("SNP",
2384
-					     is.snp,
2385
-					     is.autosome,
2386
-					     is.annotated,
2387
-					     is.X)
2388
-		marker.list <- splitIndicesByLength(marker.index, ocProbesets())
2389
-		emit <- ocLapply(seq_along(marker.list),
2390
-				 posteriorMean.snp,
2391
-				 object=object,
2392
-				 index.list=marker.list,
2393
-				 verbose=verbose,
2394
-				 prior.prob=prior.prob,
2395
-				 CN=CN,
2396
-				 scale.sd=scale.sd)
2397
-		#for(i in seq_along(marker.list)){
2398
-		#index <- marker.list[[i]]
2399
-
2400
-	#}
2401
-	} else stop("type not available")
2402
-	if(length(emit)==1) emit <- emit[[1]] else {
2403
-		state.names <- dimnames(emit[[1]])[[3]]
2404
-		tmp <- array(NA, dim=c(length(unlist(marker.list)), ncol(object), length(prior.prob)))
2405
-		for(i in seq_along(marker.list)){
2406
-			marker.index <- marker.list[[i]]
2407
-			tmp[marker.index, , ] <- emit[[i]]
2408
-		}
2409
-		emit <- tmp
2410
-		dimnames(emit) <- list(featureNames(object)[unlist(marker.list)],
2411
-				       sampleNames(object),
2412
-				       state.names)
2413
-	}#stop("need to rbind elements of emit list?")
2414
-	##tmp <- do.call("rbind", emit)
2415
-	match.index <- match(rownames(emit), featureNames(object))
2416
-	S <- length(prior.prob)
2417
-	pM <- matrix(0, length(match.index), ncol(object), dimnames=list(featureNames(object)[match.index], sampleNames(object)))
2418
-	for(i in 1:S) pM <- pM + CN[i]*emit[, , i]
2419
-	posteriorMean(object)[match.index, ] <- pM
2420
-	return(object)
2421
-}
2349
+##calculatePosteriorMean <- function(object, type=c("SNP", "NP", "X.SNP", "X.NP"), verbose=TRUE,
2350
+##				   prior.prob=c(1/7, 1/7, 3/7, 1/7, 1/7),
2351
+##				   CN=0:4, scale.sd=1){
2352
+##	stopifnot(type %in% c("SNP", "NP", "X.SNP", "X.NP"))
2353
+##	stopifnot(sum(prior.prob)==1)
2354
+##	stopifnot(length(CN) == length(prior.prob))
2355
+##	batch <- batch(object)
2356
+##	is.snp <- isSnp(object)
2357
+##	is.autosome <- chromosome(object) < 23
2358
+##	is.annotated <- !is.na(chromosome(object))
2359
+##	is.X <- chromosome(object) == 23
2360
+##	is.lds <- is(calls(object), "ffdf") | is(calls(object), "ff_matrix")
2361
+##	if(is.lds) stopifnot(isPackageLoaded("ff"))
2362
+##	samplesPerBatch <- table(as.character(batch(object)))
2363
+##	if(!"posteriorMean" %in% assayDataElementNames(object)){
2364
+##		message("adding <posteriorMean> slot to assayData")
2365
+##		pM <- matrix(NA, nrow(object), ncol(object), dimnames=list(featureNames(object), sampleNames(object)))
2366
+##		tmp <- assayDataNew(alleleA=A(object),
2367
+##				    alleleB=B(object),
2368
+##				    call=calls(object),
2369
+##				    callProbability=snpCallProbability(object),
2370
+##				    posteriorMean=pM)
2371
+##		assayData(object) <- tmp
2372
+##	}
2373
+##	## add new assay data element for posterior probabilities
2374
+##	mylabel <- function(marker.type){
2375
+##		switch(marker.type,
2376
+##		       SNP="autosomal SNPs",
2377
+##		       NP="autosomal nonpolymorphic markers",
2378
+##		       X.SNP="chromosome X SNPs",
2379
+##		       X.NP="chromosome X nonpolymorphic markers")
2380
+##	}
2381
+##	if(type=="SNP"){
2382
+##		if(verbose) message(paste("...", mylabel(type)))
2383
+##		marker.index <- whichMarkers("SNP",
2384
+##					     is.snp,
2385
+##					     is.autosome,
2386
+##					     is.annotated,
2387
+##					     is.X)
2388
+##		marker.list <- splitIndicesByLength(marker.index, ocProbesets())
2389
+##		emit <- ocLapply(seq_along(marker.list),
2390
+##				 posteriorMean.snp,
2391
+##				 object=object,
2392
+##				 index.list=marker.list,
2393
+##				 verbose=verbose,
2394
+##				 prior.prob=prior.prob,
2395
+##				 CN=CN,
2396
+##				 scale.sd=scale.sd)
2397
+##		#for(i in seq_along(marker.list)){
2398
+##		#index <- marker.list[[i]]
2399
+##
2400
+##	#}
2401
+##	} else stop("type not available")
2402
+##	if(length(emit)==1) emit <- emit[[1]] else {
2403
+##		state.names <- dimnames(emit[[1]])[[3]]
2404
+##		tmp <- array(NA, dim=c(length(unlist(marker.list)), ncol(object), length(prior.prob)))
2405
+##		for(i in seq_along(marker.list)){
2406
+##			marker.index <- marker.list[[i]]
2407
+##			tmp[marker.index, , ] <- emit[[i]]
2408
+##		}
2409
+##		emit <- tmp
2410
+##		dimnames(emit) <- list(featureNames(object)[unlist(marker.list)],
2411
+##				       sampleNames(object),
2412
+##				       state.names)
2413
+##	}#stop("need to rbind elements of emit list?")
2414
+##	##tmp <- do.call("rbind", emit)
2415
+##	match.index <- match(rownames(emit), featureNames(object))
2416
+##	S <- length(prior.prob)
2417
+##	pM <- matrix(0, length(match.index), ncol(object), dimnames=list(featureNames(object)[match.index], sampleNames(object)))
2418
+##	for(i in 1:S) pM <- pM + CN[i]*emit[, , i]
2419
+##	posteriorMean(object)[match.index, ] <- pM
2420
+##	return(object)
2421
+##}
2422 2422
 
2423 2423
 .open <- function(object){
2424 2424
 		open(B(object))
... ...
@@ -2529,6 +2529,29 @@ numberGenotypes <- function(CT){
2529 2529
 	return(res)
2530 2530
 }
2531 2531
 
2532
+
2533
+
2534
+
2535
+
2536
+genotypeCorrelation <- function(genotype){
2537
+	switch(genotype,
2538
+	       NULL="NULL",
2539
+	       A="AA",
2540
+	       B="BB",
2541
+	       AA="AA",
2542
+	       AB="AB",
2543
+	       BB="BB",
2544
+	       AAA="AA",
2545
+	       AAB="AB",
2546
+	       ABB="AB",
2547
+	       BBB="BB",
2548
+	       AAAA="AA",
2549
+	       AAAB="AB",
2550
+	       AABB="AB",
2551
+	       ABBB="AB",
2552
+	       BBBB="BB")
2553
+}
2554
+
2532 2555
 posteriorMean.snp <- function(stratum, object, index.list, CN,
2533 2556
 			      prior.prob=c(1/7, 1/7, 3/7, 1/7, 1/7), is.lds=TRUE, verbose=TRUE,
2534 2557
 			      scale.sd=1){
... ...
@@ -2651,3 +2674,28 @@ posteriorMean.snp <- function(stratum, object, index.list, CN,
2651 2674
 ##  res2[["SKW"]] <- res[["SKW"]]
2652 2675
 ##  return(list2SnpSet(res2, returnParams=returnParams))
2653 2676
 ##}
2677
+
2678
+genotypes <- function(copyNumber){
2679
+	stopifnot(copyNumber %in% 0:4)
2680
+	cn <- paste("x", copyNumber, sep="")
2681
+	switch(cn,
2682
+	       x0="NULL",
2683
+	       x1=LETTERS[1:2],
2684
+	       x2=c("AA", "AB", "BB"),
2685
+	       x3=c("AAA", "AAB", "ABB", "BBB"),
2686
+	       x4=c("AAAA", "AAAB", "AABB", "ABBB", "BBBB"))
2687
+}
2688
+
2689
+dbvn <- function(x, mu, Sigma){
2690
+	##stopifnot(ncol(x)==2)
2691
+	logA <- x[, , "A"]
2692
+	logB <- x[, , "B"]
2693
+	sigma2A <- Sigma[,1]
2694
+	sigma2B <- Sigma[,3]
2695
+	sigmaAB <- sqrt(sigma2A*sigma2B)
2696
+	rho <- Sigma[,2]
2697
+	muA <- mu[, "A"]
2698
+	muB <- mu[, "B"]
2699
+	Q.x.y <- 1/(1-rho^2)*((logA - muA)^2/sigma2A + (logB - muB)^2/sigma2B - (2*rho*((logA - muA)*(logB - muB)))/sigmaAB)
2700
+	f.x.y <- 1/(2*pi*sigmaAB*sqrt(1-rho^2))*exp(-0.5*Q.x.y)
2701
+}
... ...
@@ -37,7 +37,8 @@ setMethod("corr", signature(object="AssayData"),
37 37
 		  if(!missing.j) j <- list(...)[["j"]]
38 38
 		  batchnames <- sampleNames(object)
39 39
 		  if(!missing.j) batchnames <- batchnames[j]
40
-		  if(missing.i & missing.j) stop("Must specify either the rows i or batches j")
40
+
41
+		  ##if(missing.i & missing.j) stop("Must specify either the rows i or batches j")
41 42
 		  is.ff <- is(assayDataElement(object, "corrAA"), "ff") | is(assayDataElement(object, "corrAA"), "ffdf")
42 43
 		  if(is.ff){
43 44
 			  open(assayDataElement(object, "corrAA"))
... ...
@@ -20,294 +20,6 @@ linearParamElementReplace <- function(obj, elt, value) {
20 20
 }
21 21
 
22 22
 ## parameters
23
-setMethod("nuA", signature=signature(object="CNSet"), function(object) nu(object, "A"))
24
-setMethod("nuB", signature=signature(object="CNSet"), function(object) nu(object, "B"))
25
-setMethod("phiA", signature=signature(object="CNSet"), function(object) phi(object, "A"))
26
-setMethod("phiB", signature=signature(object="CNSet"), function(object) phi(object, "B"))
27
-setMethod("phiPrimeA", signature=signature(object="CNSet"), function(object) {
28
-	assayDataElement(batchStatistics(object), "phiPrimeA")
29
-})
30
-setMethod("phiPrimeB", signature=signature(object="CNSet"), function(object) {
31
-	assayDataElement(batchStatistics(object), "phiPrimeB")
32
-})
33
-setMethod("tau2A", signature=signature(object="CNSet"), function(object) tau2(object, "A"))
34
-setMethod("tau2B", signature=signature(object="CNSet"), function(object) tau2(object, "B"))
35
-
36
-
37
-setMethod("Ns", signature=signature(object="CNSet"),
38
-	   function(object, ...){
39
-		   Ns(batchStatistics(object), ...)
40
-	   })
41
-setMethod("corr", signature=signature(object="CNSet"),
42
-	   function(object, ...){
43
-		   corr(batchStatistics(object), ...)
44
-	   })
45
-setMethod("medians", signature=signature(object="CNSet"),
46
-	   function(object, ...){
47
-		   medians(batchStatistics(object), ...)
48
-	   })
49
-setMethod("mads", signature=signature(object="CNSet"),
50
-	   function(object, ...){
51
-		   mads(batchStatistics(object), ...)
52
-	   })
53
-setMethod("tau2", signature=signature(object="CNSet"),
54
-	   function(object, ...){
55
-		   tau2(batchStatistics(object), ...)
56
-	   })
57
-##---------------------------------------------------------------------------
58
-## Number of samples with Genotype AA, AB, or BB by batch
59
-setMethod("N.AA", signature=signature(object="CNSet"), function(object){
60
-	assayDataElement(batchStatistics(object), "N.AA")
61
-})
62
-setMethod("N.AB", signature=signature(object="CNSet"), function(object){
63
-	assayDataElement(batchStatistics(object), "N.AB")
64
-})
65
-setMethod("N.BB", signature=signature(object="CNSet"), function(object){
66
-	assayDataElement(batchStatistics(object), "N.BB")
67
-})
68
-
69
-setReplaceMethod("N.AA", signature=signature(object="CNSet", value="ff_or_matrix"),
70
-	  function(object, value){
71
-		  linearParamElementReplace(object, "N.AA", value)
72
-	  })
73
-setReplaceMethod("N.AB", signature=signature(object="CNSet", value="ff_or_matrix"),
74
-	  function(object, value){
75
-		  linearParamElementReplace(object, "N.AB", value)
76
-	  })
77
-setReplaceMethod("N.BB", signature=signature(object="CNSet", value="ff_or_matrix"),
78
-	  function(object, value){
79
-		  linearParamElementReplace(object, "N.BB", value)
80
-	  })
81
-
82
-
83
-##---------------------------------------------------------------------------
84
-##  median intensity by genotype cluster for each allele
85
-##should we update the entire matrix rather than one column...
86
-setMethod("medianA.AA", signature=signature(object="CNSet"), function(object) {
87
-	assayDataElement(batchStatistics(object), "medianA.AA")
88
-})
89
-setMethod("medianA.AB", signature=signature(object="CNSet"), function(object) {
90
-	assayDataElement(batchStatistics(object), "medianA.AB")
91
-})
92
-setMethod("medianA.BB", signature=signature(object="CNSet"), function(object) {
93
-	assayDataElement(batchStatistics(object), "medianA.BB")
94
-})
95
-setMethod("medianB.AA", signature=signature(object="CNSet"), function(object) {
96
-	assayDataElement(batchStatistics(object), "medianB.AA")
97
-})
98
-setMethod("medianB.AB", signature=signature(object="CNSet"), function(object) {
99
-	assayDataElement(batchStatistics(object), "medianB.AB")
100
-})
101
-setMethod("medianB.BB", signature=signature(object="CNSet"), function(object) {
102
-	assayDataElement(batchStatistics(object), "medianB.BB")
103
-})
104
-setReplaceMethod("medianA.AA", signature=signature(object="CNSet", value="ff_or_matrix"),
105
-	  function(object, value){
106
-		  linearParamElementReplace(object, "medianA.AA", value)
107
-	  })
108
-setReplaceMethod("medianA.AB", signature=signature(object="CNSet", value="ff_or_matrix"),
109
-	  function(object, value){
110
-		  linearParamElementReplace(object, "medianA.AB", value)
111
-	  })
112
-setReplaceMethod("medianA.BB", signature=signature(object="CNSet", value="ff_or_matrix"),
113
-	  function(object, value){
114
-		  linearParamElementReplace(object, "medianA.BB", value)
115
-	  })
116
-setReplaceMethod("medianB.AA", signature=signature(object="CNSet", value="ff_or_matrix"),
117
-	  function(object, value){
118
-		  linearParamElementReplace(object, "medianB.AA", value)
119
-	  })
120
-setReplaceMethod("medianB.AB", signature=signature(object="CNSet", value="ff_or_matrix"),
121
-	  function(object, value){
122
-		  linearParamElementReplace(object, "medianB.AB", value)
123
-	  })
124
-setReplaceMethod("medianB.BB", signature=signature(object="CNSet", value="ff_or_matrix"),
125
-	  function(object, value){
126
-		  linearParamElementReplace(object, "medianB.BB", value)
127
-	  })
128
-
129
-##---------------------------------------------------------------------------
130
-##  mad intensity by genotype cluster for each allele
131
-setMethod("madA.AA", signature=signature(object="CNSet"), function(object) {
132
-	assayDataElement(batchStatistics(object), "madA.AA")
133
-})
134
-setMethod("madA.AB", signature=signature(object="CNSet"), function(object) {
135
-	assayDataElement(batchStatistics(object), "madA.AB")
136
-})
137
-setMethod("madA.BB", signature=signature(object="CNSet"), function(object) {
138
-	assayDataElement(batchStatistics(object), "madA.BB")
139
-})
140
-setMethod("madB.AA", signature=signature(object="CNSet"), function(object) {
141
-	assayDataElement(batchStatistics(object), "madB.AA")
142
-})
143
-setMethod("madB.AB", signature=signature(object="CNSet"), function(object) {
144
-	assayDataElement(batchStatistics(object), "madB.AB")
145
-})
146
-setMethod("madB.BB", signature=signature(object="CNSet"), function(object) {
147
-	assayDataElement(batchStatistics(object), "madB.BB")
148
-})
149
-setReplaceMethod("madA.AA", signature=signature(object="CNSet", value="ff_or_matrix"),
150
-	  function(object, value){
151
-		  linearParamElementReplace(object, "madA.AA", value)
152
-	  })
153
-setReplaceMethod("madA.AB", signature=signature(object="CNSet", value="ff_or_matrix"),
154
-	  function(object, value){
155
-		  linearParamElementReplace(object, "madA.AB", value)
156
-	  })
157
-setReplaceMethod("madA.BB", signature=signature(object="CNSet", value="ff_or_matrix"),
158
-	  function(object, value){
159
-		  linearParamElementReplace(object, "madA.BB", value)
160
-	  })
161
-setReplaceMethod("madB.AA", signature=signature(object="CNSet", value="ff_or_matrix"),
162
-	  function(object, value){
163
-		  linearParamElementReplace(object, "madB.AA", value)
164
-	  })
165
-setReplaceMethod("madB.AB", signature=signature(object="CNSet", value="ff_or_matrix"),
166
-	  function(object, value){
167
-		  linearParamElementReplace(object, "madB.AB", value)
168
-	  })
169
-setReplaceMethod("madB.BB", signature=signature(object="CNSet", value="ff_or_matrix"),
170
-	  function(object, value){
171
-		  linearParamElementReplace(object, "madB.BB", value)
172
-	  })
173
-
174
-##---------------------------------------------------------------------------
175
-##  mad of log(intensities) by genotype cluster for each allele
176
-setMethod("tau2A.AA", signature=signature(object="CNSet"), function(object) {
177
-	assayDataElement(batchStatistics(object), "tau2A.AA")
178
-})
179
-setMethod("tau2A.BB", signature=signature(object="CNSet"), function(object) {
180
-	assayDataElement(batchStatistics(object), "tau2A.BB")
181
-})
182
-setMethod("tau2B.AA", signature=signature(object="CNSet"), function(object) {
183
-	assayDataElement(batchStatistics(object), "tau2B.AA")
184
-})
185
-setMethod("tau2B.BB", signature=signature(object="CNSet"), function(object) {
186
-	assayDataElement(batchStatistics(object), "tau2B.BB")
187
-})
188
-setReplaceMethod("tau2A.AA", signature=signature(object="CNSet", value="ff_or_matrix"),
189
-	  function(object, value){
190
-		  linearParamElementReplace(object, "tau2A.AA", value)
191
-	  })
192
-setReplaceMethod("tau2A.BB", signature=signature(object="CNSet", value="ff_or_matrix"),
193
-	  function(object, value){
194
-		  linearParamElementReplace(object, "tau2A.BB", value)
195
-	  })
196
-setReplaceMethod("tau2B.AA", signature=signature(object="CNSet", value="ff_or_matrix"),
197
-	  function(object, value){
198
-		  linearParamElementReplace(object, "tau2B.AA", value)
199
-	  })
200
-setReplaceMethod("tau2B.BB", signature=signature(object="CNSet", value="ff_or_matrix"),
201
-	  function(object, value){
202
-		  linearParamElementReplace(object, "tau2B.BB", value)
203
-	  })
204
-
205
-## correlation of log2A and log2B within each genotype cluster
206
-setMethod("corrAA", signature=signature(object="CNSet"), function(object){
207
-	assayDataElement(batchStatistics(object), "corrAA")
208
-  })
209
-setMethod("corrAB", signature=signature(object="CNSet"), function(object){
210
-	assayDataElement(batchStatistics(object), "corrAB")
211
-  })
212
-setMethod("corrBB", signature=signature(object="CNSet"), function(object){
213
-	assayDataElement(batchStatistics(object), "corrBB")
214
-  })
215
-setReplaceMethod("corrAA", signature=signature(object="CNSet", value="ff_or_matrix"),
216
-	  function(object, value){
217
-		  linearParamElementReplace(object, "corrAA", value)
218
-})
219
-
220
-setReplaceMethod("corrAB", signature=signature(object="CNSet", value="ff_or_matrix"),
221
-	  function(object, value){
222
-		  linearParamElementReplace(object, "corrAB", value)
223
-})
224
-
225
-setReplaceMethod("corrBB", signature=signature(object="CNSet", value="ff_or_matrix"),
226
-	  function(object, value){
227
-		  linearParamElementReplace(object, "corrBB", value)
228
-})
229
-
230
-setReplaceMethod("phiPrimeA", signature=signature(object="CNSet", value="ff_or_matrix"),
231
-	  function(object, value){
232
-		  linearParamElementReplace(object, "phiPrimeA", value)
233
-})
234
-
235
-setReplaceMethod("phiPrimeB", signature=signature(object="CNSet", value="ff_or_matrix"),
236
-	  function(object, value){
237
-		  linearParamElementReplace(object, "phiPrimeB", value)
238
-})
239
-
240
-##setMethod("mad.AA", signature=signature(object="CNSet"), function(object) mads(object, "AA"))
241
-##setMethod("mad.AB", signature=signature(object="CNSet"), function(object) mads(object, "AB"))
242
-##setMethod("mad.BB", signature=signature(object="CNSet"), function(object) mads(object, "BB"))
243
-##
244
-
245
-##setReplaceMethod("median.AA", signature=signature(object="CNSet", value="ff_or_matrix"),
246
-##	  function(object, value){
247
-##		  linearParamElementReplace(object, "median.AA", value)
248
-##	  })
249
-##setReplaceMethod("median.AB", signature=signature(object="CNSet", value="ff_or_matrix"),
250
-##	  function(object, value){
251
-##		  linearParamElementReplace(object, "median.AB", value)
252
-##	  })
253
-##setReplaceMethod("median.BB", signature=signature(object="CNSet", value="ff_or_matrix"),
254
-##	  function(object, value){
255
-##		  linearParamElementReplace(object, "median.BB", value)
256
-##	  })
257
-##setReplaceMethod("mad.AA", signature=signature(object="CNSet", value="ff_or_matrix"),
258
-##	  function(object, value){
259
-##		  linearParamElementReplace(object, "mad.AA", value)
260
-##	  })
261
-##setReplaceMethod("mad.AB", signature=signature(object="CNSet", value="ff_or_matrix"),
262
-##	  function(object, value){
263
-##		  linearParamElementReplace(object, "mad.AB", value)
264
-##	  })
265
-##setReplaceMethod("mad.BB", signature=signature(object="CNSet", value="ff_or_matrix"),
266
-##	  function(object, value){
267
-##		  linearParamElementReplace(object, "mad.BB", value)
268
-##	  })
269
-
270
-setReplaceMethod("nuA", signature=signature(object="CNSet", value="ff_or_matrix"),
271
-	  function(object, value){
272
-		  linearParamElementReplace(object, "nuA", value)
273
-	  })
274
-
275
-setReplaceMethod("nuB", signature=signature(object="CNSet", value="ff_or_matrix"),
276
-	  function(object, value){
277
-		  linearParamElementReplace(object, "nuB", value)
278
-})
279
-
280
-setReplaceMethod("phiA", signature=signature(object="CNSet", value="ff_or_matrix"),
281
-	  function(object, value){
282
-		  linearParamElementReplace(object, "phiA", value)
283
-})
284
-
285
-setReplaceMethod("phiB", signature=signature(object="CNSet", value="ff_or_matrix"),
286
-	  function(object, value){
287
-		  linearParamElementReplace(object, "phiB", value)
288
-})
289
-
290
-setReplaceMethod("tau2A", signature=signature(object="CNSet", value="ff_or_matrix"),
291
-	  function(object, value){
292
-		  linearParamElementReplace(object, "tau2A", value)
293
-})
294
-
295
-setReplaceMethod("tau2B", signature=signature(object="CNSet", value="ff_or_matrix"),
296
-	  function(object, value){
297
-		  linearParamElementReplace(object, "tau2B", value)
298
-})
299
-
300
-
301
-
302
-setReplaceMethod("flags", signature=signature(object="CNSet", value="ff_or_matrix"),
303
-		 function(object, value){
304
-			 linearParamElementReplace(object, "flags", value)
305
-})
306
-setReplaceMethod("flags", signature=signature(object="CNSet", value="ff_matrix"),
307
-		 function(object, value){
308
-			 linearParamElementReplace(object, "flags", value)
309
-})
310
-
311 23
 
312 24
 ## allele A
313 25
 ##   autosome SNPs
... ...
@@ -571,22 +283,101 @@ setMethod("totalCopynumber", signature=signature(object="CNSet"),
571 283
 	  })
572 284
 rawCopynumber <- totalCopynumber
573 285
 
574
-##setReplaceMethod("snpCall", c("CNSet", "ff_or_matrix"),
575
-##                 function(object, ..., value){
576
-##			 assayDataElementReplace(object, "call", value)
577
-##		 })
286
+setMethod("posteriorProbability", signature(object="CNSet"),
287
+	  function(object, predictRegion, copyNumber=0:4){
288
+		  logI <- array(NA, dim=c(nrow(object), ncol(object), 2), dimnames=list(NULL, NULL, LETTERS[1:2]))
289
+		  getIntensity <- function(object){
290
+			  logI[, , 1] <- log2(A(object))
291
+			  logI[, , 2] <- log2(B(object))
292
+			  return(logI)
293
+		  }
294
+		  logI <- getIntensity(object)
295
+		  ##gts <- lapply(as.list(0:4), genotypes)
296
+		  prob <- array(NA, dim=c(nrow(object), ncol(object), length(copyNumber)),
297
+				dimnames=list(NULL,NULL, paste("copynumber",copyNumber, sep="")))
298
+		  for(i in seq_along(copyNumber)){
299
+			  G <- genotypes(copyNumber[i])
300
+			  P <- array(NA, dim=c(nrow(object), ncol(object), length(G)),
301
+				     dimnames=list(NULL,NULL,G))
302
+			  for(g in seq_along(G)){
303
+				  gt <- G[g]
304
+				  P[, , gt] <- dbvn(x=logI, mu=predictRegion[[gt]]$mu, Sigma=predictRegion[[gt]]$cov)
305
+			  }
306
+			  ## the marginal probability for the total copy number
307
+			  ##  -- integrate out the probability for the different genotypes
308
+			  for(j in 1:ncol(P)){
309
+				  prob[, j, i] <- rowSums(as.matrix(P[, j, ]), na.rm=TRUE)
310
+			  }
311
+		  }
312
+		  ## divide by normalizing constant.  Probs across copy number states must sum to one.
313
+		  ##nc <- apply(prob, c(1,3), sum, na.rm=TRUE)
314
+		  nc <- matrix(NA, nrow(object), ncol(object))
315
+		  for(j in seq_len(ncol(object))){
316
+			  nc <- rowSums(as.matrix(prob[,j, ]), na.rm=TRUE)
317
+			  prob[, j, ] <- prob[, j, ]/nc
318
+		  }
319
+		  return(prob)
320
+	  })
578 321
 
579
-setReplaceMethod("snpCall", "CNSet",
580
-                 function(object, value){
581
-			 assayDataElementReplace(object, "call", value)
582
-		 })
322
+setMethod("calculatePosteriorMean", signature(object="CNSet"),
323
+	  function(object, posteriorProb, copyNumber=0:4, w, ...){
324
+		  if(missing(w)) w <- rep(1/length(copyNumber),length(copyNumber))
325
+		  stopifnot(sum(w)==1)
326
+		  stopifnot(dim(posteriorProb)[[3]] == length(copyNumber))
327
+		  pm <- matrix(0, nrow(object), ncol(object))
328
+		  for(i in seq_along(copyNumber)){
329
+			  pm <- pm + posteriorProb[, , i] * copyNumber[i] * w[i]
330
+		  }
331
+		  return(pm)
332
+	  })
583 333
 
584
-##setReplaceMethod("snpCallProbability", c("CNSet", "ff_or_matrix"),
585
-##                 function(object, ..., value){
586
-##			 assayDataElementReplace(object, "callProbability", value)
587
-##		 })
588 334
 
589
-setReplaceMethod("snpCallProbability", "CNSet",
590
-                 function(object, value){
591
-			 assayDataElementReplace(object, "callProbability", value)
592
-		 })
335
+## for a given copy number, return a named list of bivariate normal prediction regions
336
+##   - elements of list are named by genotype
337
+##   'AAA'  list should have 'mu' and 'Sigma'
338
+##                mu is a R x 2 matrix, R is number of features
339
+##                Sigma is a R x 4 matrix.  Each row can be coerced to a 2 x 2 matrix
340
+##          For nonpolymorphic markers, 2nd column of mu is NA and elements 2-4 of Sigma are NA
341
+##
342
+setMethod("predictionRegion", signature(object="CNSet", copyNumber="integer"),
343
+	  function(object, copyNumber=0:4){
344
+		  ## would it be more efficient to store as an array?
345
+		  ## mu: features x genotype x allele
346
+		  ## Sigma: features x genotype x covariance
347
+		  stopifnot(all(copyNumber %in% 0:4))
348
+		  gts <- lapply(as.list(copyNumber), genotypes)
349
+		  nms <- unlist(gts)
350
+		  res <- vector("list", length(nms))
351
+		  ##names(res) <- paste("copyNumber", copyNumber, sep="")
352
+		  names(res) <- nms
353
+		  nu <- getNu(object)
354
+		  phi <- getPhi(object)
355
+		  mus <- matrix(NA, nrow(nu), 2, dimnames=list(NULL, LETTERS[1:2]))
356
+		  Sigma <- matrix(NA, nrow(nu), 3)
357
+		  bivariateCenter <- function(nu, phi){
358
+			  ##  lexical scope for mus, CA, CB
359
+			  mus[,1] <- log2(nu[1] + CA * phi[1])
360
+			  mus[,2] <- log2(nu[2] + CB * phi[2])
361
+			  mus
362
+		  }
363
+		  for(i in seq_along(copyNumber)){
364
+			  G <- genotypes(copyNumber[i])
365
+			  tmp <- vector("list", length(G))
366
+			  names(tmp) <- G
367
+			  CN <- copyNumber[i]
368
+			  for(g in seq_along(G)){
369
+				  gt <- G[g]
370
+				  CB <- g-1
371
+				  CA <- CN-CB
372
+				  gt.corr <- genotypeCorrelation(gt)
373
+				  nma <- ifelse(CA == 0, "tau2A.BB", "tau2A.AA")
374
+				  nmb <- ifelse(CB == 0, "tau2B.AA", "tau2B.BB")
375
+				  Sigma[,1] <- getVar(object, nma)
376
+				  Sigma[,3] <- getVar(object, nmb)
377
+				  Sigma[,2] <- getCor(object, gt.corr)
378
+				  res[[gt]]$mu <- bivariateCenter(nu, phi)
379
+				  res[[gt]]$cov <- Sigma
380
+			  }
381
+		  }
382
+		  return(res)
383
+	  })