Browse code

added a few tests for totalCopyNumber in the illumina_copynumber vignette

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

Rob Scharp authored on 30/08/2010 15:45:37
Showing4 changed files

... ...
@@ -72,3 +72,4 @@ export(crlmm,
72 72
 ##export(initializeParamObject, biasAdjNP2)
73 73
 
74 74
 
75
+exportMethods(totalCopyNumber)
75 76
\ No newline at end of file
... ...
@@ -157,12 +157,12 @@ ellipse.CNSet <- function(x, copynumber, batch, ...){
157 157
 
158 158
 
159 159
 setMethod("totalCopyNumber",
160
-	  signature=signature(object="CNSet", ...),
160
+	  signature=signature(object="CNSet"),
161 161
 	  function(object, ...){
162 162
 		  is.ff <- is(CA(object), "ff") | is(CA(object), "ffdf")
163
-		  dotArgs <- names(list(...))
164
-		  missing.i <- "i" %in% dotArgs
165
-		  missing.j <- "j" %in% dotArgs
163
+		  dotArgs <- list(...)
164
+		  missing.i <- !"i" %in% names(dotArgs)
165
+		  missing.j <- !"j" %in% names(dotArgs)
166 166
 		  if(missing.i & missing.j){
167 167
 			  if(is.ff) stop("Must specify i and/or j for ff objects")
168 168
 		  }
... ...
@@ -173,7 +173,8 @@ setMethod("totalCopyNumber",
173 173
 			  snp.index2 <- which(isSnp(object)[i])
174 174
 		  } else {
175 175
 			  i <- 1:nrow(object)
176
-			  snp.index2 <- snp.index <- which(isSnp(object))
176
+			  snp.index <- which(isSnp(object))
177
+			  snp.index2 <- snp.index
177 178
 		  }
178 179
 		  if(!missing.j){
179 180
 			  j <- dotArgs[["j"]]
... ...
@@ -181,8 +182,7 @@ setMethod("totalCopyNumber",
181 182
 		  cn.total <- as.matrix(CA(object)[i, j])
182 183
 		  if(length(snp.index) > 0){
183 184
 			  cb <- as.matrix(CB(object)[snp.index, j])
184
-			  ##				 snps <- (1:nrow(cn.total))[i %in% snp.index]
185
-			  cn.total[snp.index2, ] <- cn.total[snp.index, ] + cb
185
+			  cn.total[snp.index2, ] <- cn.total[snp.index2, ] + cb
186 186
 		  }
187 187
 		  cn.total <- cn.total/100
188 188
 		  return(cn.total)
... ...
@@ -319,6 +319,8 @@ genotypeConf <- integerScoreToProbability(snpCallProbability(x)[snp.index[1:10],
319 319
 Allele-specific copy number at polymorphic loci:
320 320
 <<ca>>=
321 321
 ca <- CA(x[snp.index, ])/100
322
+cb <- CB(x[snp.index, ])/100
323
+ct <- ca+cb
322 324
 @
323 325
 
324 326
 Total copy number at nonpolymorphic loci:
... ...
@@ -328,9 +330,9 @@ cn.nonpolymorphic <- CA(x[np.index, ])/100
328 330
 
329 331
 Total copy number at both polymorphic and nonpolymorphic loci:
330 332
 <<totalCopynumber>>=
331
-path <- system.file("scripts", package="crlmm")
332
-source(file.path(path, "helperFunctions.R"))
333
-cn <- totalCopyNumber(x, i=c(snp.index, np.index))
333
+##path <- system.file("scripts", package="crlmm")
334
+##source(file.path(path, "helperFunctions.R"))
335
+cn <- crlmm:::totalCopyNumber(x, i=c(snp.index, np.index))
334 336
 apply(cn, 2, median, na.rm=TRUE)
335 337
 @
336 338
 
... ...
@@ -251,23 +251,27 @@ cb <- cb[-missing.index, ]
251 251
 \noindent Negating the \Robject{isSnp} function could be used to
252 252
 extract the estimates at nonpolymorphic markers. For instance,
253 253
 <<monomorphic-accessor>>=
254
-cn.monomorphic <- CA(cnSet)[which(chromosome(cnSet) == 21 & !isSnp(cnSet)), ]/100
254
+np.index <- which(chromosome(cnSet) == 21 & !isSnp(cnSet))
255
+cn.monomorphic <- CA(cnSet)[np.index, ]/100
255 256
 @
256 257
 
257 258
 At polymorphic loci, the total copy number is the sum of the number of
258 259
 copies of the A allele and the number of copies for the B allele.  At
259 260
 nonpolymorphic loci, the total copy number is the number of copies for
260
-the A allele.  The helper function totalCopyNumber can be used to
261
-extract the total copy number for all polymorphic and nonpolymorphic
262
-markers.  Documentation of the \Rfunction{totalCopyNumber} will be
263
-available in the devel version of the \Rpackage{oligoClasses}.
261
+the A allele.  The helper function \Rfunction{totalCopyNumber} can be
262
+used to extract the total copy number for all polymorphic and
263
+nonpolymorphic markers.  Documentation of the
264
+\Rfunction{totalCopyNumber} will be available in the devel version of
265
+the \Rpackage{oligoClasses}.
264 266
 
265
-<<copyNumberHelper>>=
267
+<<totalCopyNumber>>=
266 268
 cn.total <- ca+cb
267
-cn.total2 <- totalCopyNumber(cnSet, i=marker.index)
269
+cn.total2 <- crlmm:::totalCopyNumber(cnSet, i=marker.index)
268 270
 cn.total2 <- cn.total2[-missing.index, ]
269
-dimnames(cn.total) <- NULL
270
-all.equal(cn.total2, cn.total)
271
+stopifnot(all.equal(cn.total2, cn.total))
272
+
273
+cn.monomorphic2 <- totalCopyNumber(cnSet, i=np.index)
274
+stopifnot(all.equal(cn.monomorphic, cn.monomorphic2))
271 275
 @
272 276
 
273 277
 A few simple visualizations may be helpful at this point. The first
... ...
@@ -322,40 +326,6 @@ axis(1, at=pretty(x), labels=pretty(x/1e6))
322 326
 mtext("Mb", 1, outer=TRUE, line=2)
323 327
 @
324 328
 
325
-%Next we remove some of the waves using \Rpackage{limma}'s function
326
-%\Rfunction{loessFit}.
327
-%<<wavecorrection, fig=TRUE>>=
328
-%require(limma)
329
-%par(mfrow=c(2,2), las=1, mar=c(0.5, 4, 0.5, 0.5), oma=c(4, 1, 1,1))
330
-%for(j in c(low.snr, high.snr)){
331
-%	y <- copyNumber(cnSet)[chromosome(cnSet) == 1, j]/100
332
-%	missing <- is.na(y)
333
-%	x <- x[!missing]
334
-%	y <- y[!missing]
335
-%	dist <- diff(x)
336
-%	index <- c(0, cumsum(log10(dist) > 6))
337
-%	x.split <- split(x, index)
338
-%	y <- as.numeric(runmed(y, k=3))
339
-%	y.split <- split(y, index)
340
-%	y.smooth <- vector("list", length(y.split))
341
-%	for(i in seq_along(y.smooth)){
342
-%		fit <- loessFit(y=y.split[[i]], x= x.split[[i]], span=0.3)
343
-%		y.smooth[[i]] <- 2+fit$residuals
344
-%	}
345
-%	plot(unlist(x.split),
346
-%	     unlist(y.split),
347
-%	     pch=".", ylab="copy number", xaxt="n", log="y", ylim=c(0.5,5))
348
-%	plot(unlist(x.split),
349
-%	     unlist(y.smooth),
350
-%	     pch=".", ylab="copy number", xaxt="n", log="y", ylim=c(0.5,5), yaxt="n")
351
-%	legend("topright", bty="n", legend=paste("SNR =", round(cnSet$SNR[j], 1)))
352
-%	abline(h=2, col="grey70")
353
-%	if(j == high.snr){
354
-%		axis(1, at=pretty(x), labels=pretty(x/1e6))
355
-%		mtext("Mb", 1, outer=TRUE, line=2)
356
-%	}
357
-%}
358
-%@
359 329
 \section{Session information}
360 330
 <<sessionInfo, results=tex>>=
361 331
 toLatex(sessionInfo())