Browse code

fix for finding index of nonpolymorphic markers in genotype function

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

Rob Scharp authored on 16/02/2011 15:56:48
Showing 3 changed files

... ...
@@ -47,7 +47,7 @@ getFeatureData <- function(cdfName, copynumber=FALSE){
47 47
 		M[index, "isSnp"] <- 0L
48 48
 	}
49 49
 	##A few of the snpProbes do not match -- I think it is chromosome Y.
50
-	M <- M[!is.na(M[, "chromosome"]), ]
50
+	M[is.na(M[, "isSnp"]), "isSnp"] <- 1L
51 51
 	M <- data.frame(M)
52 52
 	return(new("AnnotatedDataFrame", data=M))
53 53
 }
... ...
@@ -163,8 +163,16 @@ if(!file.exists(file.path(outdir, "cnSet.rda"))){
163 163
 			     .FUN=genotype,
164 164
 			     filenames=celFiles,
165 165
 			     cdfName=cdfName,
166
-			     batch=batch, .load.it=FALSE)
167
-	class(calls(gtSet))
166
+			     batch=batch)
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 176
 }
169 177
 @
170 178
 
... ...
@@ -189,8 +197,7 @@ The \Rfunction{crlmmCopynumber} performs the following steps:
189 197
   displayed during the processing.
190 198
 
191 199
 <<LDS_copynumber>>=
192
-GT.CONF.THR <- 0.8
193
-cnSet <- checkExists("cnSet", .path=outdir, .FUN=crlmmCopynumber, object=gtSet, GT.CONF.THR=GT.CONF.THR)
200
+cnSet <- checkExists("cnSet", .path=outdir, .FUN=crlmmCopynumber, object=gtSet, .load.it=FALSE)
194 201
 @
195 202
 
196 203
 In an effort to reduce I/O, the \Rpackage{crlmmCopynumber} function no
... ...
@@ -219,7 +226,7 @@ the useful accessors for extracting which markers are SNPs, which are
219 226
 chromosomes, etc.
220 227
 
221 228
 <<ca>>=
222
-snp.index <- which(isSnp(cnSet) & chromosome(cnSet) < 23)
229
+snp.index <- which(isSnp(cnSet) & !is.na(chromosome(cnSet)))
223 230
 ca <- CA(cnSet, i=snp.index, j=1:5)
224 231
 cb <- CB(cnSet, i=snp.index, j=1:5)
225 232
 ct <- ca+cb
... ...
@@ -244,20 +251,22 @@ ct2 <- totalCopynumber(cnSet, i=marker.index, j=1:5)
244 251
 stopifnot(all.equal(ct, ct2))
245 252
 @
246 253
 
247
-  Nonpolymorphic markers on chromosome X:
248 254
 
249
-<<nonpolymorphicX, fig=TRUE, width=8, height=4>>=
250
-library(RColorBrewer)
251
-bp.cols <- brewer.pal(8, "Paired")[3:4]
255
+TODO: FIX estimation for nonpolymorphic markers on X
256
+
257
+<<nonpolymorphicX>>=
252 258
 npx.index <- which(chromosome(cnSet)==23 & !isSnp(cnSet))
259
+a <- A(cnSet)[npx.index, M]
260
+nus <- nuA(cnSet)[npx.index, ]
261
+## nu and phi are not estimated
253 262
 M <- sample(which(cnSet$gender==1), 5)
254
-cols <- bp.cols[cnSet$gender[c(M,F)]]
255
-cn.M <- CA(cnSet, i=npx.index, j=M)
256
-cn.F <- CA(cnSet, i=npx.index, j=F)
257 263
 F <- sample(which(cnSet$gender==2), 5)
258
-par(las=1)
259
-boxplot(data.frame(cbind(cn.M, cn.F)), pch=".", col=cols, outline=FALSE)
260
-legend("topleft", bty="n", fill=bp.cols, legend=c("male", "female"))
264
+X.markers <- which(!isSnp(cnSet) & chromosome(cnSet) == 23)
265
+cn.M <- CA(cnSet, i=X.markers, j=M)
266
+cn.F <- CA(cnSet, i=X.markers, j=F)
267
+##phi(cnSet, "A")[X.markers[1:10]]
268
+index <- which(is.na(nuA(cnSet)[X.markers, 1]))
269
+boxplot(data.frame(cbind(cn.M, cn.F)), pch=".")
261 270
 @
262 271
 
263 272
 \begin{figure}[t!]
... ...
@@ -1,3 +1,2 @@
1
-rsync -avuzb --exclude '*~' --exclude '.git*' -e ssh ~/madman/Rpacks/git/crlmm/ enigma2.jhsph.edu:~/madman/Rpacks/git/crlmm
2
-rsync -avuzb --exclude '*~' --exclude '.git*' -e ssh enigma2.jhsph.edu:~/madman/Rpacks/git/crlmm/inst/scripts/*opynumber.pdf ~/madman/Rpacks/git/crlmm/inst/scripts/
1
+rsync -avuzb --exclude '.git*' -e ssh ~/madman/Rpacks/git/crlmm/ enigma2.jhsph.edu:~/madman/Rpacks/release/crlmm
3 2