Browse code

Updates to genotype function and fit.lm4

genotype function uses ffcolapply and assigns results from snprma to
specific indices. This change was needed so that genotype() is not
dependendent on a specific version of the genomewidesnp*Crlmm
annotation package.

Changes to fit.lm4 fix estimates of copy number at nonpolymorphic loci
on chromosome X.

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

Rob Scharp authored on 15/11/2010 14:21:05
Showing1 changed files

... ...
@@ -46,6 +46,7 @@ getFeatureData <- function(cdfName, copynumber=FALSE){
46 46
 		M[index, "chromosome"] <- chromosome2integer(cnProbes[, grep("chr", colnames(cnProbes))])
47 47
 		M[index, "isSnp"] <- 0L
48 48
 	}
49
+	##A few of the snpProbes do not match -- I think it is chromosome Y.
49 50
 	M <- M[!is.na(M[, "chromosome"]), ]
50 51
 	M <- data.frame(M)
51 52
 	return(new("AnnotatedDataFrame", data=M))
... ...
@@ -145,6 +146,13 @@ genotype <- function(filenames,
145 146
 	##message("Saving snprmaRes file")
146 147
 	##save(snprmaRes, file=file.path(outdir, "snprmaRes.rda"))
147 148
 	if(verbose) message("Finished preprocessing.")
149
+	gns <- snprmaRes[["gns"]]
150
+	snp.I <- isSnp(callSet)
151
+	is.snp <- which(snp.I)
152
+	snp.index <- match(featureNames(callSet)[is.snp], gns)
153
+	stopifnot(identical(featureNames(callSet)[is.snp], gns[snp.index]))
154
+	if(is.lds) open(callSet)
155
+	mixtureParams <- matrix(NA, 4, length(filenames))
148 156
 	if(is.lds){
149 157
 		open(snprmaRes[["A"]])
150 158
 		open(snprmaRes[["B"]])
... ...
@@ -155,24 +163,15 @@ genotype <- function(filenames,
155 163
 			   BATCHBYTES=bb)
156 164
 		ffcolapply(B(callSet)[is.snp, i1:i2] <- snprmaRes[["B"]][snp.index, i1:i2], X=snprmaRes[["B"]],
157 165
 			   BATCHBYTES=bb)
158
-##		##bb <- getOption("ffbatchbytes")
159
-##		message("Writing normalized intensities to callSet")
160
-##		system.time(
161
-##		for(j in 1:ncol(callSet)){
162
-##			A(callSet)[is.snp, j] <- snprmaRes[["A"]][snp.index, j]
163
-##			B(callSet)[is.snp, j] <- snprmaRes[["B"]][snp.index, j]
164
-##		})
165
-		##ffrowapply(A(callSet)[i1:i2, ] <- snprmaRes[["A"]][i1:i2, ], X=snprmaRes[["A"]])##, BATCHBYTES=bb)
166
-		##ffrowapply(B(callSet)[i1:i2, ] <- snprmaRes[["B"]][i1:i2, ], X=snprmaRes[["B"]])##, BATCHBYTES=bb)
167
-		pData(callSet)$SKW <- snprmaRes[["SKW"]]
168
-		pData(callSet)$SNR <- snprmaRes[["SNR"]]
169 166
 	} else{
170 167
 		A(callSet)[is.snp, ] <- snprmaRes[["A"]][snp.index, ]
171 168
 		B(callSet)[is.snp, ] <- snprmaRes[["B"]][snp.index, ]
172
-		pData(callSet)$SKW <- snprmaRes[["SKW"]]
173
-		pData(callSet)$SNR <- snprmaRes[["SNR"]]
174 169
 	}
175
-	np.index <- which(!snp.I)
170
+	pData(callSet)$SKW <- snprmaRes[["SKW"]]
171
+	pData(callSet)$SNR <- snprmaRes[["SNR"]]
172
+	mixtureParams <- snprmaRes$mixtureParams
173
+	np.index <- which(!is.snp)
174
+	if(verbose) message("Normalizing nonpolymorphic markers")
176 175
 	FUN <- ifelse(is.lds, "cnrma2", "cnrma")
177 176
 	## main purpose is to update 'alleleA'
178 177
 	cnrmaFxn <- function(FUN,...){
... ...
@@ -222,18 +221,20 @@ genotype <- function(filenames,
222 221
 			   BATCHBYTES=bb)
223 222
 		ffcolapply(snpCallProbability(callSet)[is.snp, i1:i2] <- tmp[["B"]][snp.index, i1:i2], X=tmp[["B"]],
224 223
 			   BATCHBYTES=bb)
225
-##		for(j in 1:ncol(callSet)){
226
-##			snpCall(callSet)[is.snp, j] <- tmp[["calls"]][snp.index, j]
227
-##			snpCallProbability(callSet)[is.snp, j] <- tmp[["confs"]][snp.index, j]
228
-##		}
229
-		##		ffrowapply(snpCall(callSet)[i1:i2, ] <- tmp[["calls"]][i1:i2, ], X=tmp[["calls"]])#, BATCHBYTES=bb)
230
-		##		ffrowapply(snpCallProbability(callSet)[i1:i2, ] <- tmp[["confs"]][i1:i2, ], X=tmp[["confs"]])#, BATCHBYTES=bb)
224
+		close(tmp[["calls"]])
225
+		close(tmp[["confs"]])
231 226
 	} else {
232 227
 		calls(callSet)[is.snp, ] <- tmp[["calls"]][snp.index, ]
233 228
 		snpCallProbability(callSet)[is.snp, ] <- tmp[["confs"]][snp.index, ]
234 229
 	}
235 230
 	message("Finished updating.  Cleaning up.")
236 231
 	callSet$gender <- tmp$gender
232
+	if(is.lds){
233
+		delete(snprmaRes[["A"]])
234
+		delete(snprmaRes[["B"]])
235
+		delete(snprmaRes[["mixtureParams"]])
236
+		rm(snprmaRes)
237
+	}
237 238
 	close(callSet)
238 239
 	if(is.lds){
239 240
 		delete(snprmaRes[["A"]])
... ...
@@ -936,13 +937,11 @@ fit.lm4 <- function(strata,
936 937
 		Y <- log2(c(tmp[, 3], tmp[, 4]))
937 938
 		betahat <- solve(crossprod(X), crossprod(X, Y))
938 939
 		if(enough.men){
939
-			##X.men <- cbind(1, medianA.AA.M[, k])
940 940
 			X.men <- cbind(1, log2(medianA.AA.M[, k]))
941 941
 			Yhat1 <- as.numeric(X.men %*% betahat)
942 942
 			## put intercept and slope for men in nuB and phiB
943 943
 			phiB[, k] <- 2^(Yhat1)
944
-			##nuB[, k] <- 2^(medianA.AA.M[, k]) - phiB[, k]
945
-			nuB[, k] <- medianA.AA.M[, k] - 1*phiB[, k] ## male has 1 copy
944
+			nuB[, k] <- medianA.AA.M[, k] - 1*phiB[, k]
946 945
 		}
947 946
 		X.fem <- cbind(1, log2(medianA.AA.F[, k]))
948 947
 		Yhat2 <- as.numeric(X.fem %*% betahat)