Browse code

added ACN method in methods-CNSet.R. CA and CB call ACN.

- need to check for chromosome X,Y

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

Rob Scharp authored on 21/08/2010 02:48:26
Showing3 changed files

... ...
@@ -76,7 +76,7 @@ export(constructIlluminaCNSet)
76 76
 export(linesCNSetLM)
77 77
 export(computeCN, fit.lm1, fit.lm2, fit.lm3, fit.lm4, construct,
78 78
        dqrlsWrapper, nuphiAllele)
79
-export(computeCopynumber)
79
+export(computeCopynumber, ACN)
80 80
 
81 81
 
82 82
 
... ...
@@ -237,36 +237,59 @@ setMethod("corr", c("CNSetLM", "character"), function(object, allele){
237 237
 	return(res)
238 238
 })
239 239
 
240
-setMethod("CB", 
241
-	  signature=signature(object="CNSet", i="integerOrMissing", j="integerOrMissing"),
242
-	  function(object, i, j){
243
-		  ##assayDataElement(object, "CB")
244
-		  browser()
245
-	  })
240
+ACN <- function(object, allele, i , j){
241
+	if(missing(i) & missing(j)){
242
+		if(inherits(A(object), "ff") | inherits(A(object), "ffdf")) stop("Must specify i and/or j for ff objects")
243
+	}
244
+	if(missing(i) & !missing(j)){
245
+		## calculate ca only for batches indexed by j
246
+		ubatch <- unique(batch(object))
247
+		batches <- unique(batch(object)[j])
248
+		for(k in seq_along(batches)){
249
+			l <- match(batches[k], ubatch)
250
+			bg <- nu(object, allele)[, l]
251
+			sl <- phi(object, allele)[, l]
252
+			I <- allele(object, allele)[, j]
253
+			acn <- 1/sl*(I - bg)				  
254
+		}
255
+	}
256
+	if(!missing(i) & missing(j)){
257
+		## calculate ca, cb for all batches
258
+		batches <- unique(batch(object))
259
+		for(k in seq_along(batches)){
260
+			##bb <- batches[k]
261
+			bg <- nu(object, allele)[i, k]
262
+			sl <- phi(object, allele)[i, k]
263
+			I <- allele(object, allele)[i, j]
264
+			acn <- 1/sl*(I - bg)
265
+		}
266
+	}
267
+	if(!missing(i) & !missing(j)){
268
+		ubatch <- unique(batch(object))
269
+		batches <- unique(batch(object)[j])
270
+		for(k in seq_along(batches)){
271
+			l <- match(batches[k], ubatch)
272
+			bg <- nu(object, allele)[i, l]
273
+			sl <- phi(object, allele)[i, l]
274
+			I <- allele(object, allele)[i, j]
275
+			acn <- 1/bg*(I - sl)				  
276
+		}			  
277
+	}
278
+	return(acn)
279
+}
246 280
 
247 281
 setMethod("CA",
248 282
 	  signature=signature(object="CNSet", i="integerOrMissing", j="integerOrMissing"),
249 283
 	  function(object, i, j) {
250
-		  browser()
251
-		  if(missing(i) & missing(j)){
252
-			  if(inherits(A(object), "ff") | inherits(A(object), "ffdf")) stop("Must specify i and/or j for ff objects")
253
-		  }
254
-		  if(missing(i) & !missing(j)){
255
-			  ## calculate ca only for batches indexed by j
256
-			  batches <- unique(batch(object))[j]
257
-			  for(k in seq_along(batches)){
258
-				  bb <- batches[k]
259
-			  }
260
-		  }
261
-		  if(!missing(i) & missing(j)){
262
-			  ## calculate ca, cb for all batches
263
-			  batches <- unique(batch(object))
264
-		  }
265
-		  if(!missing(i) & !missing(j)){
266
-
267
-		  }
284
+		  ca <- ACN(object, allele="A", i, j)
268 285
 		  return(ca)
269 286
 	  })
287
+setMethod("CB",
288
+	  signature=signature(object="CNSet", i="integerOrMissing", j="integerOrMissing"),
289
+	  function(object, i, j) {
290
+		  cb <- ACN(object, allele="B", i, j)
291
+		  return(cb)
292
+	  })
270 293
 
271 294
 setMethod("totalCopyNumber",
272 295
 	  signature=signature(object="CNSet", i="integerOrMissing", j="integerOrMissing"),
... ...
@@ -301,7 +324,7 @@ setMethod("totalCopyNumber",
301 324
 			cn.total[snps, ] <- cn.total[snps, ] + cb
302 325
 		}
303 326
 	}
304
-	cn.total <- cn.total/100
327
+	##cn.total <- cn.total/100
305 328
 	dimnames(cn.total) <- NULL
306 329
 	return(cn.total)
307 330
 })
... ...
@@ -130,7 +130,6 @@ if(!file.exists(file.path(outdir, "cnSet.rda"))){
130 130
 	class(calls(gtSet.assayData_matrix))
131 131
 }
132 132
 ##q("no")
133
-stop()
134 133
 @ 
135 134
 
136 135
 Next, we estimate copy number for the 20 CEL files.  The copy number
... ...
@@ -147,9 +146,9 @@ cnSet.assayData_matrix <- checkExists("cnSet.assayData_matrix",
147 146
 				      .FUN=crlmmCopynumber,
148 147
 				      object=gtSet.assayData_matrix,
149 148
 				      chromosome=22)
150
-Rprof(interval=0.1)
151
-obj <- crlmmCopynumber(gtSet.assayData_matrix, chromosome=22)
152
-Rprof(NULL)
149
+##Rprof(interval=0.1)
150
+##obj <- crlmmCopynumber(gtSet.assayData_matrix, chromosome=22)
151
+##Rprof(NULL)
153 152
 if(file.exists(file.path(outdir, "gtSet.assayData_matrix.rda")))
154 153
 	unlink(file.path(outdir, "gtSet.assayData_matrix.rda"))
155 154
 @ 
... ...
@@ -181,13 +180,14 @@ so that subsequent calls to \verb@Sweave@ can be run interactively.
181 180
 
182 181
 <<LDS_genotype>>=
183 182
 if(!file.exists(file.path(outdir, "cnSet.assayData_ff.rda"))){
184
-	library(ff)
183
+	Rprof(filename="Rprof_genotypeff.out", interval=0.1)
185 184
 	gtSet.assayData_ff <- checkExists("gtSet.assayData_ff",
186 185
 					  .path=outdir,
187 186
 					  .FUN=genotypeLD,
188 187
 					  filenames=celFiles,
189 188
 					  cdfName=cdfName,
190 189
 					  batch=batch)
190
+	Rprof(NULL)
191 191
 	class(calls(gtSet.assayData_ff))
192 192
 }
193 193
 @ 
... ...
@@ -198,12 +198,14 @@ practice, once the object \Robject{cnSet.assayData_ff} is created the
198 198
 from the \Robject{outdir}.
199 199
 
200 200
 <<LDS_copynumber>>=
201
+Rprof(filename="Rprof_cnff.out", interval=0.1)
201 202
 cnSet.assayData_ff <- checkExists("cnSet.assayData_ff",
202 203
 				  .path=outdir,
203 204
 				  .FUN=crlmmCopynumberLD,
204 205
 				  filenames=celFiles,
205 206
 				  cdfName=cdfName,
206 207
 				  batch=batch)
208
+Rprof(NULL)
207 209
 if(file.exists(file.path(outdir, "gtSet.assayData_ff.rda")))
208 210
 	unlink(file.path(outdir, "gtSet.assayData_ff.rda"))
209 211
 gc()
... ...
@@ -373,7 +375,11 @@ genotypeConf <- integerScoreToProbability(snpCallProbability(x)[snp.index[1:10],
373 375
 
374 376
 Allele-specific copy number at polymorphic loci:
375 377
 <<ca>>=
376
-ca <- CA(x[snp.index, ])/100
378
+##ca <- CA(x[snp.index, ])/100
379
+snp.index <- which(isSnp(obj))
380
+ca <- CA(obj, i=snp.index)
381
+##or
382
+ca <- ACN(obj, "A", i=snp.index)
377 383
 @ 
378 384
 
379 385
 Total copy number at nonpolymorphic loci: