Browse code

Fix bug in CB accessor for nonpolymorphic markers

CB returns 0 at nonpolymorphic loci instead of NA. While NA makes more
sense since we do not measure the alternative allele, this prevents
NA's from creeping into estimates of total copy number (CA+CB = NA vs
CA+CB=CA)

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

Rob Scharp authored on 15/11/2010 14:21:52
Showing 2 changed files

... ...
@@ -326,7 +326,7 @@ C1 <- function(object, marker.index, batch.index, sample.index){
326 326
 	return(as.matrix(acn))
327 327
 }
328 328
 
329
-## allele B  (allele 'A' for chromosome X NPs)
329
+## allele B  (treated allele 'A' for chromosome X NPs)
330 330
 ##   autosome SNPs
331 331
 ##   chromosome X for male nonpolymorphic markers
332 332
 C2 <- function(object, marker.index, batch.index, sample.index, NP.X=FALSE){
... ...
@@ -493,9 +493,11 @@ ACN <- function(object, allele, i , j){
493 493
 			open(B(object))
494 494
 		}
495 495
 		if(any(is.auto)){
496
-			auto.index <- which(is.auto)
497
-			marker.index <- i[is.auto]
498
-			acn[auto.index, ] <- C2(object, marker.index, batch.index, j)
496
+			auto.index <- which(is.auto & is.snp)
497
+			if(length(auto.index) > 0){
498
+				marker.index <- i[is.auto]
499
+				acn[auto.index, ] <- C2(object, marker.index, batch.index, j)
500
+			}
499 501
 		}
500 502
 		if(any(is.X)){
501 503
 			##Chr X SNPs
... ...
@@ -165,14 +165,11 @@ if(!file.exists(file.path(outdir, "cnSet.rda"))){
165 165
 			     cdfName=cdfName,
166 166
 			     batch=batch)
167 167
 	## try normalizing np probes on chr X
168
-	callSet <- gtSet
169
-	filenames <- celFiles
170
-	snp.I <- isSnp(callSet)
171
-	is.snp <- which(snp.I)
172
-	np.index <- which(!is.snp)
173
-
174
-
175
-
168
+##	callSet <- gtSet
169
+##	filenames <- celFiles
170
+##	snp.I <- isSnp(callSet)
171
+##	is.snp <- which(snp.I)
172
+##	np.index <- which(!is.snp)
176 173
 }
177 174
 @
178 175
 
... ...
@@ -253,22 +250,20 @@ ct2 <- totalCopynumber(cnSet, i=marker.index, j=1:5)
253 250
 stopifnot(all.equal(ct, ct2))
254 251
 @
255 252
 
256
-Nonpolymorphic markers on X:
257 253
 
258
-\begin{figure}[!t]
259
-\centering
260
-<<nonpolymorphicX>>=
254
+Nonpolymorphic markers on chromosome X:
255
+
256
+\begin{figure}
257
+  \centering
258
+<<nonpolymorphicX, fig=TRUE>>=
261 259
 npx.index <- which(chromosome(cnSet)==23 & !isSnp(cnSet))
262
-a <- A(cnSet)[npx.index, M]
263
-nus <- nuA(cnSet)[npx.index, ]
264
-## nu and phi are not estimated
265 260
 M <- sample(which(cnSet$gender==1), 5)
266 261
 F <- sample(which(cnSet$gender==2), 5)
267
-X.markers <- which(!isSnp(cnSet) & chromosome(cnSet) == 23)
268
-cn.M <- CA(cnSet, i=X.markers, j=M)
269
-cn.F <- CA(cnSet, i=X.markers, j=F)
270
-##phi(cnSet, "A")[X.markers[1:10]]
271
-index <- which(is.na(nuA(cnSet)[X.markers, 1]))
262
+##a <- A(cnSet)[npx.index, M]
263
+##nus <- nuA(cnSet)[npx.index, ]
264
+## nu and phi are not estimated
265
+cn.M <- CA(cnSet, i=npx.index, j=M)
266
+cn.F <- CA(cnSet, i=npx.index, j=F)
272 267
 boxplot(data.frame(cbind(cn.M, cn.F)), pch=".")
273 268
 @
274 269
 \caption{Copy number estimates for nonpolymorphic markers on
... ...
@@ -277,13 +272,13 @@ boxplot(data.frame(cbind(cn.M, cn.F)), pch=".")
277 272
   males and 2 for females.}
278 273
 \end{figure}
279 274
 
280
-
281 275
 Polymorphic markers, X chromosome:
282 276
 
283 277
 \begin{figure}
284 278
   \centering
285 279
 <<polymorphicX, fig=TRUE, width=8, height=4>>=
286
-X.markers <- sample(which(isSnp(cnSet) & chromosome(cnSet) == 23), 25e3)
280
+## copy number estimates on X for SNPs are biased towards small values.
281
+X.markers <- which(isSnp(cnSet) & chromosome(cnSet) == 23)
287 282
 ca.M <- CA(cnSet, i=X.markers, j=M)
288 283
 cb.M <- CB(cnSet, i=X.markers, j=M)
289 284
 ca.F <- CA(cnSet, i=X.markers, j=F)
... ...
@@ -294,8 +289,9 @@ boxplot(data.frame(cbind(cn.M, cn.F)), pch=".", outline=FALSE, col=cols, xaxt="n
294 289
 legend("topleft", fill=unique(cols), legend=c("Male", "Female"))
295 290
 abline(h=c(1,2))
296 291
 ## alternatively
297
-cn.F2 <- totalCopynumber(cnSet, i=X.markers, j=F)
298
-stopifnot(all.equal(cn.F, cn.F2))
292
+cn.MF1 <- cbind(cn.M, cn.F)
293
+cn.MF2 <- totalCopynumber(cnSet, i=X.markers, j=c(M,F))
294
+##stopifnot(all.equal(cn.MF1[!is.na(cn.MF1)], cn.MF2[!is.na(cn.MF2)]))
299 295
 @
300 296
 \caption{Copy number estimates for polymorphic markers on chromosome
301 297
   X.  We assume that the median copy number (across samples) at any