Browse code

methods Ns, medians, corr, and mads return arrays.

Ns and corr return an array with dimension: features x genotype x batch

medians and mads return an array with dimension: features x allele x genotype x batch

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

Rob Scharp authored on 21/08/2010 02:50:19
Showing5 changed files

... ...
@@ -76,7 +76,7 @@ export(crlmm,
76 76
 export(constructIlluminaCNSet)
77 77
 export(totalCopynumber)
78 78
 export(cnrma, cnrma2)
79
-exportMethods(A, B, nuA, nuB, phiA, phiB, corr, tau2, Ns)
79
+exportMethods(A, B, nuA, nuB, phiA, phiB, corr, tau2, Ns, medians, mads)
80 80
 export(genotypeSummary,
81 81
        shrinkSummary,
82 82
        estimateCnParameters)
... ...
@@ -1,99 +1,198 @@
1 1
 setMethod("Ns", signature(object="AssayData"),
2
-	  function(object, genotype, batchname){
3
-		  if(missing(batchname) & length(sampleNames(object)) > 1)
4
-			  stop("must supply batchname")
5
-		  if(!missing(batchname)){
6
-			  stopifnot(batchname %in% batchNames(object))
7
-			  j <- match(batchname, sampleNames(object))
8
-		  } else j <- 1
9
-		  if(missing(genotype)){
10
-			  res <- cbind(assayDataElement(object, "N.AA")[, j],
11
-				       assayDataElement(object, "N.AB")[, j],
12
-				       assayDataElement(object, "N.BB")[, j])
13
-			  return(res)
14
-		  } else{
15
-			  getValue <- function(genotype){
16
-				  switch(genotype,
17
-					 AA="N.AA",
18
-					 AB="N.AB",
19
-					 BB="N.BB",
20
-					 stop("allele must be 'AA', 'AB', or 'BB'"))
21
-			  }
22
-			  val <- getValue(genotype)
23
-			  return(assayDataElement(object, val)[, j])
2
+	  function(object, i, j, ...){
3
+		  if(!missing(j)){
4
+			  batchnames <- unique(as.character(batch(object)[j]))
5
+		  } else batchnames <- batchNames(object)
6
+		  nc <- length(batchnames)
7
+		  if(!missing(i)) nr <- length(i) else nr <- nrow(object)
8
+		  res <- array(NA, dim=c(nr, 3, nc))
9
+		  dimnames(res)[[2]] <- c("AA", "AB", "BB")
10
+		  dimnames(res)[[3]] <- batchnames
11
+		  if(missing(i) & missing(j)){
12
+			  N.AA <- as.matrix(assayDataElement(object, "N.AA"))
13
+			  N.AB <- as.matrix(assayDataElement(object, "N.AB"))
14
+			  N.BB <- as.matrix(assayDataElement(object, "N.BB"))
15
+		  }
16
+		  if(missing(i) & !missing(j)){
17
+			  J <- match(batchnames, batchNames(object))
18
+			  stopifnot(length(J) > 0 & !all(is.na(J)))
19
+			  N.AA <- as.matrix(assayDataElement(object, "N.AA"))[, J, ...]
20
+			  N.AB <- as.matrix(assayDataElement(object, "N.AB"))[, J, ...]
21
+			  N.BB <- as.matrix(assayDataElement(object, "N.BB"))[, J, ...]
24 22
 		  }
23
+		  if(!missing(i) & !missing(j)){
24
+			  J <- match(batchnames, batchNames(object))
25
+			  stopifnot(length(J) > 0 & !all(is.na(J)))
26
+			  N.AA <- as.matrix(assayDataElement(object, "N.AA"))[i, J, ...]
27
+			  N.AB <- as.matrix(assayDataElement(object, "N.AB"))[i, J, ...]
28
+			  N.BB <- as.matrix(assayDataElement(object, "N.BB"))[i, J, ...]
29
+		  }
30
+		  if(!missing(i) & missing(j)){
31
+			  N.AA <- as.matrix(assayDataElement(object, "N.AA"))[i, , ...]
32
+			  N.AB <- as.matrix(assayDataElement(object, "N.AB"))[i, , ...]
33
+			  N.BB <- as.matrix(assayDataElement(object, "N.BB"))[i, , ...]
34
+		  }
35
+		  res[, "AA", ] <- N.AA
36
+		  res[, "AB", ] <- N.AB
37
+		  res[, "BB", ] <- N.BB
38
+		  return(res)
25 39
 	  })
26 40
 setMethod("corr", signature(object="AssayData"),
27
-	  function(object, genotype, batchname){
28
-		  if(missing(batchname) & length(sampleNames(object)) > 1)
29
-			  stop("must supply batchname")
30
-		  if(!missing(batchname)){
31
-			  stopifnot(batchname %in% batchNames(object))
32
-			  j <- match(batchname, sampleNames(object))
33
-		  } else j <- 1
34
-		  if(missing(genotype)){
35
-			  res <- cbind(assayDataElement(object, "corrAA")[, j],
36
-				       assayDataElement(object, "corrAB")[, j],
37
-				       assayDataElement(object, "corrBB")[, j])
38
-			  return(res)
39
-		  } else{
40
-			  getValue <- function(genotype){
41
-				  switch(genotype,
42
-					 AA="corrAA",
43
-					 AB="corrAB",
44
-					 BB="corrBB",
45
-					 stop("allele must be 'AA', 'AB', or 'BB'"))
46
-			  }
47
-			  val <- getValue(genotype)
48
-			  return(assayDataElement(object, val)[, j])
41
+	  function(object, i, j, ...){
42
+		  if(!missing(j)){
43
+			  batchnames <- unique(as.character(batch(object)[j]))
44
+		  } else batchnames <- batchNames(object)
45
+		  nc <- length(batchnames)
46
+		  if(!missing(i)) nr <- length(i) else nr <- nrow(object)
47
+		  res <- array(NA, dim=c(nr, 3, nc))
48
+		  dimnames(res)[[2]] <- c("AA", "AB", "BB")
49
+		  dimnames(res)[[3]] <- batchnames
50
+		  if(missing(i) & missing(j)){
51
+			  corrAA <- as.matrix(assayDataElement(object, "corrAA"))
52
+			  corrAB <- as.matrix(assayDataElement(object, "corrAB"))
53
+			  corrBB <- as.matrix(assayDataElement(object, "corrBB"))
54
+		  }
55
+		  if(missing(i) & !missing(j)){
56
+			  J <- match(batchnames, batchNames(object))
57
+			  stopifnot(length(J) > 0 & !all(is.na(J)))
58
+			  corrAA <- as.matrix(assayDataElement(object, "corrAA"))[, J, ...]
59
+			  corrAB <- as.matrix(assayDataElement(object, "corrAB"))[, J, ...]
60
+			  corrBB <- as.matrix(assayDataElement(object, "corrBB"))[, J, ...]
61
+		  }
62
+		  if(!missing(i) & !missing(j)){
63
+			  J <- match(batchnames, batchNames(object))
64
+			  stopifnot(length(J) > 0 & !all(is.na(J)))
65
+			  corrAA <- as.matrix(assayDataElement(object, "corrAA"))[i, J, ...]
66
+			  corrAB <- as.matrix(assayDataElement(object, "corrAB"))[i, J, ...]
67
+			  corrBB <- as.matrix(assayDataElement(object, "corrBB"))[i, J, ...]
49 68
 		  }
69
+		  if(!missing(i) & missing(j)){
70
+			  corrAA <- as.matrix(assayDataElement(object, "corrAA"))[i, , ...]
71
+			  corrAB <- as.matrix(assayDataElement(object, "corrAB"))[i, , ...]
72
+			  corrBB <- as.matrix(assayDataElement(object, "corrBB"))[i, , ...]
73
+		  }
74
+		  res[, "AA", ] <- corrAA
75
+		  res[, "AB", ] <- corrAB
76
+		  res[, "BB", ] <- corrBB
77
+		  return(res)
50 78
 	  })
51 79
 
52 80
 setMethod("medians", signature(object="AssayData"),
53
-	  function(object, allele, batchname){
54
-		  stopifnot(!missing(allele))
55
-		  if(missing(batchname) & length(sampleNames(object)) > 1)
56
-			  stop("must supply batchname")
57
-		  if(!missing(batchname)){
58
-			  stopifnot(batchname %in% batchNames(object))
59
-			  j <- match(batchname, sampleNames(object))
60
-		  } else j <- 1
61
-		  getMedians <- function(allele){
62
-			  switch(allele,
63
-				 A=cbind(assayDataElement(object, "medianA.AA")[, j],
64
-				         assayDataElement(object, "medianA.AB")[, j],
65
-         			         assayDataElement(object, "medianA.BB")[, j]),
66
-				 B=cbind(assayDataElement(object, "medianB.AA")[, j],
67
-				         assayDataElement(object, "medianB.AB")[, j],
68
-				         assayDataElement(object, "medianB.BB")[, j]),
69
-				 stop("allele must be 'A' or 'B'")
70
-				 )
81
+	  function(object, i, j, ...){
82
+		  if(!missing(j)){
83
+			  batchnames <- unique(as.character(batch(object)[j]))
84
+		  } else batchnames <- batchNames(object)
85
+		  nc <- length(batchnames)
86
+		  if(!missing(i)) nr <- length(i) else nr <- nrow(object)
87
+		  res <- array(NA, dim=c(nr, 2, 3, nc))
88
+		  dimnames(res)[[2]] <- c("A", "B")
89
+		  dimnames(res)[[3]] <- c("AA", "AB", "BB")
90
+		  dimnames(res)[[4]] <- batchnames
91
+		  if(missing(i) & missing(j)){
92
+			  medianA.AA <- as.matrix(assayDataElement(object, "medianA.AA"))
93
+			  medianA.AB <- as.matrix(assayDataElement(object, "medianA.AB"))
94
+			  medianA.BB <- as.matrix(assayDataElement(object, "medianA.BB"))
95
+			  medianB.AA <- as.matrix(assayDataElement(object, "medianB.AA"))
96
+			  medianB.AB <- as.matrix(assayDataElement(object, "medianB.AB"))
97
+			  medianB.BB <- as.matrix(assayDataElement(object, "medianB.BB"))
98
+		  }
99
+		  if(missing(i) & !missing(j)){
100
+			  J <- match(batchnames, batchNames(object))
101
+			  stopifnot(length(J) > 0 & !all(is.na(J)))
102
+			  medianA.AA <- as.matrix(assayDataElement(object, "medianA.AA"))[, J, ...]
103
+			  medianA.AB <- as.matrix(assayDataElement(object, "medianA.AB"))[, J, ...]
104
+			  medianA.BB <- as.matrix(assayDataElement(object, "medianA.BB"))[, J, ...]
105
+			  medianB.AA <- as.matrix(assayDataElement(object, "medianB.AA"))[, J, ...]
106
+			  medianB.AB <- as.matrix(assayDataElement(object, "medianB.AB"))[, J, ...]
107
+			  medianB.BB <- as.matrix(assayDataElement(object, "medianB.BB"))[, J, ...]
108
+
109
+		  }
110
+		  if(!missing(i) & !missing(j)){
111
+			  J <- match(batchnames, batchNames(object))
112
+			  stopifnot(length(J) > 0 & !all(is.na(J)))
113
+			  medianA.AA <- as.matrix(assayDataElement(object, "medianA.AA"))[i, J, ...]
114
+			  medianA.AB <- as.matrix(assayDataElement(object, "medianA.AB"))[i, J, ...]
115
+			  medianA.BB <- as.matrix(assayDataElement(object, "medianA.BB"))[i, J, ...]
116
+			  medianB.AA <- as.matrix(assayDataElement(object, "medianB.AA"))[i, J, ...]
117
+			  medianB.AB <- as.matrix(assayDataElement(object, "medianB.AB"))[i, J, ...]
118
+			  medianB.BB <- as.matrix(assayDataElement(object, "medianB.BB"))[i, J, ...]
119
+		  }
120
+		  if(!missing(i) & missing(j)){
121
+			  medianA.AA <- as.matrix(assayDataElement(object, "medianA.AA"))[i, ...]
122
+			  medianA.AB <- as.matrix(assayDataElement(object, "medianA.AB"))[i, ...]
123
+			  medianA.BB <- as.matrix(assayDataElement(object, "medianA.BB"))[i, ...]
124
+			  medianB.AA <- as.matrix(assayDataElement(object, "medianB.AA"))[i, ...]
125
+			  medianB.AB <- as.matrix(assayDataElement(object, "medianB.AB"))[i, ...]
126
+			  medianB.BB <- as.matrix(assayDataElement(object, "medianB.BB"))[i, ...]
71 127
 		  }
72
-		  getMedians(allele)
128
+		  res[, "A", "AA", ] <- medianA.AA
129
+		  res[, "A", "AB", ] <- medianA.AB
130
+		  res[, "A", "BB", ] <- medianA.BB
131
+		  res[, "B", "AA", ] <- medianB.AA
132
+		  res[, "B", "AB", ] <- medianB.AB
133
+		  res[, "B", "BB", ] <- medianB.BB
134
+		  return(res)
73 135
 })
74
-setMethod("mads", signature(object="AssayData"),
75
-	  function(object, allele, batchname){
76
-		  stopifnot(!missing(allele))
77
-		  if(missing(batchname) & length(sampleNames(object)) > 1)
78
-			  stop("must supply batchname")
79
-		  if(!missing(batchname)){
80
-			  stopifnot(batchname %in% batchNames(object))
81
-			  j <- match(batchname, sampleNames(object))
82
-		  } else j <- 1
83
-		  getMads <- function(allele){
84
-			  switch(allele,
85
-				 A=cbind(assayDataElement(object, "madA.AA")[, j],
86
-				         assayDataElement(object, "madA.AB")[, j],
87
-         			         assayDataElement(object, "madA.BB")[, j]),
88
-				 B=cbind(assayDataElement(object, "madB.AA")[, j],
89
-				         assayDataElement(object, "madB.AB")[, j],
90
-				         assayDataElement(object, "madB.BB")[, j]),
91
-				 stop("allele must be 'A' or 'B'")
92
-				 )
136
+
137
+setMethod("medians", signature(object="AssayData"),
138
+	  function(object, i, j, ...){
139
+		  if(!missing(j)){
140
+			  batchnames <- unique(as.character(batch(object)[j]))
141
+		  } else batchnames <- batchNames(object)
142
+		  nc <- length(batchnames)
143
+		  if(!missing(i)) nr <- length(i) else nr <- nrow(object)
144
+		  res <- array(NA, dim=c(nr, 2, 3, nc))
145
+		  dimnames(res)[[2]] <- c("A", "B")
146
+		  dimnames(res)[[3]] <- c("AA", "AB", "BB")
147
+		  dimnames(res)[[4]] <- batchnames
148
+		  if(missing(i) & missing(j)){
149
+			  madA.AA <- as.matrix(assayDataElement(object, "madA.AA"))
150
+			  madA.AB <- as.matrix(assayDataElement(object, "madA.AB"))
151
+			  madA.BB <- as.matrix(assayDataElement(object, "madA.BB"))
152
+			  madB.AA <- as.matrix(assayDataElement(object, "madB.AA"))
153
+			  madB.AB <- as.matrix(assayDataElement(object, "madB.AB"))
154
+			  madB.BB <- as.matrix(assayDataElement(object, "madB.BB"))
93 155
 		  }
94
-		  getMads(allele)
156
+		  if(missing(i) & !missing(j)){
157
+			  J <- match(batchnames, batchNames(object))
158
+			  stopifnot(length(J) > 0 & !all(is.na(J)))
159
+			  madA.AA <- as.matrix(assayDataElement(object, "madA.AA"))[, J, ...]
160
+			  madA.AB <- as.matrix(assayDataElement(object, "madA.AB"))[, J, ...]
161
+			  madA.BB <- as.matrix(assayDataElement(object, "madA.BB"))[, J, ...]
162
+			  madB.AA <- as.matrix(assayDataElement(object, "madB.AA"))[, J, ...]
163
+			  madB.AB <- as.matrix(assayDataElement(object, "madB.AB"))[, J, ...]
164
+			  madB.BB <- as.matrix(assayDataElement(object, "madB.BB"))[, J, ...]
165
+
166
+		  }
167
+		  if(!missing(i) & !missing(j)){
168
+			  J <- match(batchnames, batchNames(object))
169
+			  stopifnot(length(J) > 0 & !all(is.na(J)))
170
+			  madA.AA <- as.matrix(assayDataElement(object, "madA.AA"))[i, J, ...]
171
+			  madA.AB <- as.matrix(assayDataElement(object, "madA.AB"))[i, J, ...]
172
+			  madA.BB <- as.matrix(assayDataElement(object, "madA.BB"))[i, J, ...]
173
+			  madB.AA <- as.matrix(assayDataElement(object, "madB.AA"))[i, J, ...]
174
+			  madB.AB <- as.matrix(assayDataElement(object, "madB.AB"))[i, J, ...]
175
+			  madB.BB <- as.matrix(assayDataElement(object, "madB.BB"))[i, J, ...]
176
+		  }
177
+		  if(!missing(i) & missing(j)){
178
+			  madA.AA <- as.matrix(assayDataElement(object, "madA.AA"))[i, ...]
179
+			  madA.AB <- as.matrix(assayDataElement(object, "madA.AB"))[i, ...]
180
+			  madA.BB <- as.matrix(assayDataElement(object, "madA.BB"))[i, ...]
181
+			  madB.AA <- as.matrix(assayDataElement(object, "madB.AA"))[i, ...]
182
+			  madB.AB <- as.matrix(assayDataElement(object, "madB.AB"))[i, ...]
183
+			  madB.BB <- as.matrix(assayDataElement(object, "madB.BB"))[i, ...]
184
+		  }
185
+		  res[, "A", "AA", ] <- madA.AA
186
+		  res[, "A", "AB", ] <- madA.AB
187
+		  res[, "A", "BB", ] <- madA.BB
188
+		  res[, "B", "AA", ] <- madB.AA
189
+		  res[, "B", "AB", ] <- madB.AB
190
+		  res[, "B", "BB", ] <- madB.BB
191
+		  return(res)
95 192
 })
96 193
 
194
+
195
+
97 196
 setMethod("tau2", signature(object="AssayData"),
98 197
 	  function(object, allele, batchname){
99 198
 		  stopifnot(!missing(allele))
... ...
@@ -421,24 +421,22 @@ setMethod("CB", signature=signature(object="CNSet"),
421 421
 		  return(cb)
422 422
 	  })
423 423
 
424
-##totalCopynumber <- function(object, ...){
425 424
 setMethod("totalCopynumber", signature=signature(object="CNSet"),
426 425
 	  function(object, ...){
427
-	message("copy number at nonpolymorphic probes is not currently available")
428
-	ca <- CA(object, ...)
429
-	cb <- CB(object, ...)
430
-	is.snp <- isSnp(object)
431
-	dotArgs <- list(...)
432
-	if("i" %in% names(dotArgs)){
433
-		i <- dotArgs[["i"]]
434
-		np.index <- which(!is.snp[i])
435
-		if(length(np.index) > 0) cb[np.index, ] <- 0
436
-	} else {
437
-		np.index <- which(!is.snp)
438
-		if(length(np.index) > 0) cb[np.index, ] <- 0
439
-	}
440
-	return(ca+cb)
441
-})
426
+		  ca <- CA(object, ...)
427
+		  cb <- CB(object, ...)
428
+		  is.snp <- isSnp(object)
429
+		  dotArgs <- list(...)
430
+		  if("i" %in% names(dotArgs)){
431
+			  i <- dotArgs[["i"]]
432
+			  np.index <- which(!is.snp[i])
433
+			  if(length(np.index) > 0) cb[np.index, ] <- 0
434
+		  } else {
435
+			  np.index <- which(!is.snp)
436
+			  if(length(np.index) > 0) cb[np.index, ] <- 0
437
+		  }
438
+		  return(ca+cb)
439
+	  })
442 440
 
443 441
 setReplaceMethod("snpCall", c("CNSet", "ff_or_matrix"),
444 442
                  function(object, ..., value){
... ...
@@ -10,7 +10,7 @@
10 10
 crlmmCopynumber(object, MIN.SAMPLES=10, SNRMin = 5, MIN.OBS = 1, 
11 11
 	        DF.PRIOR = 50, bias.adj = FALSE,
12 12
                 prior.prob = rep(1/4, 4), seed = 1, verbose = TRUE,
13
-                GT.CONF.THR = 0.95, PHI.THR = 2^6, MIN.NU = 2^3, MIN.PHI = 2^3,
13
+                GT.CONF.THR = 0.95, MIN.NU = 2^3, MIN.PHI = 2^3,
14 14
                 THR.NU.PHI = TRUE, type=c("SNP", "NP", "X.SNP", "X.NP"))
15 15
 }
16 16
 
... ...
@@ -36,7 +36,7 @@ library(ff)
36 36
 path <- system.file("extdata", package="crlmm")
37 37
 ldPath(path)
38 38
 data(sample.CNSetLMff)
39
-cnSet <- as(sample.CNSetLMff, "CNSet")
39
+cnSetff <- as(sample.CNSetLMff, "CNSet")
40 40
 all(isCurrent(cnSet))
41 41
 
42 42
 ## a bigger object with multiple batches
... ...
@@ -91,18 +91,12 @@ nu(cnSet, "A")[1:5, ]
91 91
 nu(cnSet, "B")[1:5, ]
92 92
 ## the slope
93 93
 phi(cnSet, "A")[1:5, ]
94
-## the variance for CN > 0 (log2-scale)
95
-sigma2(cnSet, "A")[1:5, ]
96
-sigma2(cnSet, "B")[1:5, ]
97
-## background variance (log2-scale)
98
-tau2(cnSet, "A")[1:5, ]
99
-tau2(cnSet, "B")[1:5, ]
100 94
 ## correlation within genotype cluster AA
101
-corr(cnSet, "AA")[1:5, ]
102
-## correlation within genotype cluster AB
103
-corr(cnSet, "AB")[1:5, ]
104
-## correlation within genotype cluster BB
105
-corr(cnSet, "BB")[1:5, ]
95
+##corr(cnSet, "AA")[1:5, ]
96
+#### correlation within genotype cluster AB
97
+##corr(cnSet, "AB")[1:5, ]
98
+#### correlation within genotype cluster BB
99
+##corr(cnSet, "BB")[1:5, ]
106 100
 ## --------------------------------------------------
107 101
 
108 102
 ## --------------------------------------------------