Browse code

fixed cyclic namespace dependency by updating version dependencies in DESCRIPTION, and updating other bioc packages

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

Rob Scharp authored on 02/03/2011 00:59:03
Showing4 changed files

... ...
@@ -9,8 +9,8 @@ Description: Faster implementation of CRLMM specific to SNP 5.0 and 6.0 arrays,
9 9
 License: Artistic-2.0
10 10
 Depends: R (>= 2.13.0),
11 11
          methods,
12
-         Biobase (>= 2.11.8),
13
-         oligoClasses (>= 1.13.9)
12
+         Biobase (>= 2.11.9),
13
+         oligoClasses (>= 1.13.11)
14 14
 Imports: affyio (>= 1.19.2),
15 15
          ellipse,
16 16
          ff (>= 2.2-1),
... ...
@@ -92,5 +92,8 @@ setGeneric("tau2A<-", function(object, value) standardGeneric("tau2A<-"))
92 92
 setGeneric("tau2B<-", function(object, value) standardGeneric("tau2B<-"))
93 93
 setGeneric("flags<-", function(object, value) standardGeneric("flags<-"))
94 94
 
95
+setGeneric("posteriorMean", function(object) standardGeneric("posteriorMean"))
96
+setGeneric("posteriorMean<-", function(object, value) standardGeneric("posteriorMean<-"))
97
+
95 98
 
96 99
 
... ...
@@ -2153,8 +2153,11 @@ imputeAcrossBatch <- function(N.AA, N.AB, N.BB,
2153 2153
 	return(list(res, updated))
2154 2154
 }
2155 2155
 
2156
-posteriorMean <- function(object, type=c("SNP", "NP", "X.SNP", "X.NP"), verbose=TRUE,
2157
-			  prior.prob=c(1/6, 1/6, 3/6, 1/6)){
2156
+
2157
+
2158
+
2159
+calculatePosteriorMean <- function(object, type=c("SNP", "NP", "X.SNP", "X.NP"), verbose=TRUE,
2160
+			  prior.prob=c(1/7, 1/7, 3/7, 1/7, 1/7)){
2158 2161
 	stopifnot(type %in% c("SNP", "NP", "X.SNP", "X.NP"))
2159 2162
 	stopifnot(sum(prior.prob)==1)
2160 2163
 	batch <- batch(object)
... ...
@@ -2165,6 +2168,8 @@ posteriorMean <- function(object, type=c("SNP", "NP", "X.SNP", "X.NP"), verbose=
2165 2168
 	is.lds <- is(calls(object), "ffdf") | is(calls(object), "ff_matrix")
2166 2169
 	if(is.lds) stopifnot(isPackageLoaded("ff"))
2167 2170
 	samplesPerBatch <- table(as.character(batch(object)))
2171
+
2172
+	## add new assay data element for posterior probabilities
2168 2173
 	mylabel <- function(marker.type){
2169 2174
 		switch(marker.type,
2170 2175
 		       SNP="autosomal SNPs",
... ...
@@ -2186,16 +2191,20 @@ posteriorMean <- function(object, type=c("SNP", "NP", "X.SNP", "X.NP"), verbose=
2186 2191
 				 index.list=marker.list,
2187 2192
 				 verbose=verbose,
2188 2193
 				 prior.prob=prior.prob)
2189
-	}
2194
+		for(i in seq_along(marker.list)){
2195
+			index <- marker.list[[i]]
2196
+
2197
+		}
2198
+	} else stop("type not available")
2190 2199
 	return(emit)
2191 2200
 }
2192 2201
 
2193 2202
 posteriorMean.snp <- function(stratum, object, index.list,
2194
-			      prior.prob=c(1/6, 1/6, 3/6, 1/6), is.lds=TRUE, verbose=TRUE){
2203
+			      prior.prob=c(1/7, 1/7, 3/7, 1/7, 1/7), is.lds=TRUE, verbose=TRUE){
2195 2204
 	if(verbose) message("Probe stratum ", stratum, " of ", length(index.list))
2196 2205
 	index <- index.list[[stratum]]
2197
-	if(is.lds){
2198
-		open(A(object))
2206
+	test <- tryCatch(open(A(object)), error=function(e) NULL)
2207
+	if(!is.null(test)){
2199 2208
 		open(B(object))
2200 2209
 		open(tau2A.AA(object))
2201 2210
 		open(tau2B.BB(object))
... ...
@@ -2211,6 +2220,7 @@ posteriorMean.snp <- function(stratum, object, index.list,
2211 2220
 	}
2212 2221
 	a <- log2(as.matrix(A(object)[index, ]))
2213 2222
 	b <- log2(as.matrix(B(object)[index, ]))
2223
+	NN <- Ns(object, i=index)[, , 1]
2214 2224
 	sig2A <- as.matrix(tau2A.AA(object)[index,])
2215 2225
 	sig2B <- as.matrix(tau2B.BB(object)[index,])
2216 2226
 	tau2A <- as.matrix(tau2A.BB(object)[index,])
... ...
@@ -2222,7 +2232,7 @@ posteriorMean.snp <- function(stratum, object, index.list,
2222 2232
 	phiA <- as.matrix(phiA(object)[index, ])
2223 2233
 	nuB <- as.matrix(nuB(object)[index, ])
2224 2234
 	phiB <- as.matrix(phiB(object)[index, ])
2225
-	if(is.lds){
2235
+	if(!is.null(test)){
2226 2236
 		close(A(object))
2227 2237
 		close(B(object))
2228 2238
 		close(tau2A.AA(object))
... ...
@@ -2237,7 +2247,8 @@ posteriorMean.snp <- function(stratum, object, index.list,
2237 2247
 		close(phiA(object))
2238 2248
 		close(phiB(object))
2239 2249
 	}
2240
-	emit <- array(NA, dim=c(nrow(a), ncol(a), 4))##SNPs x sample x 'truth'
2250
+	S <- length(prior.prob)
2251
+	emit <- array(NA, dim=c(nrow(a), ncol(a), S))##SNPs x sample x 'truth'
2241 2252
 	sample.index <- split(1:ncol(object), batch(object))
2242 2253
 	sum.mymatrix <- function(...){
2243 2254
 		x <- list(...)
... ...
@@ -2247,9 +2258,9 @@ posteriorMean.snp <- function(stratum, object, index.list,
2247 2258
 	for(j in seq_along(sample.index)){
2248 2259
 		cat("batch ", j, "\n")
2249 2260
 		J <- sample.index[[j]]
2250
-		probs <- array(NA, dim=c(nrow(a), length(J), 4))
2251
-		for(k in seq_along(0:3)){
2252
-			CT <- (0:3)[k]
2261
+		probs <- array(NA, dim=c(nrow(a), length(J), S))
2262
+		for(k in seq_along(0:4)){
2263
+			CT <- (0:4)[k]
2253 2264
 			f.x.y <- vector("list", choose(CT+1, CT))
2254 2265
 			for(i in seq_along(0:CT)){
2255 2266
 				CA <- (0:CT)[i]
... ...
@@ -2276,9 +2287,9 @@ posteriorMean.snp <- function(stratum, object, index.list,
2276 2287
 		}
2277 2288
 		emit[, J, ] <- probs
2278 2289
 	}
2279
-	for(i in 1:4) emit[, , i] <- emit[, , i]*prior.prob[i]
2290
+	for(i in 1:S) emit[, , i] <- emit[, , i]*prior.prob[i]
2280 2291
 	total <- matrix(0, nrow(emit), ncol(emit))
2281
-	for(i in 1:4) total <- total+emit[, , i]
2292
+	for(i in 1:S) total <- total+emit[, , i]
2282 2293
 	message(" need to check for outliers before rescaling")
2283 2294
 	is.outlier <- total < 1e-5
2284 2295
 ##	plot(a[outlier.index[1], ], b[outlier.index[1], ], col="grey50",
... ...
@@ -2292,7 +2303,9 @@ posteriorMean.snp <- function(stratum, object, index.list,
2292 2303
 	## how to handle outliers...
2293 2304
 	##  - use priors (posterior mean will likely be near 2).  smoothing needs to take into account the uncertainty
2294 2305
 	##  - need uncertainty estimates for posterior means
2295
-	for(i in 1:4) emit[, , i] <- emit[, , i]/total
2306
+	for(i in 1:S) emit[, , i] <- emit[, , i]/total
2307
+	dimnames(emit) <- featureNames(object)[index]
2308
+	## for one marker/one sample, the emission probs must sum to 1
2296 2309
 	return(emit)
2297 2310
 }
2298 2311
 
... ...
@@ -1,3 +1,6 @@
1
+setMethod("posteriorMean", signature(object="CNSet"), function(object) assayDataElement(object, "posteriorMean"))
2
+setReplaceMethod("posteriorMean", signature(object="CNSet", value="matrix"), function(object, value) assayDataElementReplace(object, "posteriorMean", value))
3
+
1 4
 linearParamElementReplace <- function(obj, elt, value) {
2 5
     storage.mode <- storageMode(batchStatistics(obj))
3 6
     switch(storage.mode,