Browse code

Added Ingo to the crlmm authorship. Moved Rafa to the last position.

Took care of some NAMESPACE issues -- was exporting functions that
are no longer in cnrma-functions.R

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

Rob Scharp authored on 21/08/2010 02:48:43
Showing4 changed files

... ...
@@ -3,7 +3,7 @@ Type: Package
3 3
 Title: Genotype Calling (CRLMM) and Copy Number Analysis tool for Affymetrix SNP 5.0 and 6.0 and Illumina arrays.
4 4
 Version: 1.7.11
5 5
 Date: 2010-07-30
6
-Author: Rafael A Irizarry, Benilton S Carvalho <carvalho@bclab.org>, Robert Scharpf <rscharpf@jhsph.edu>, Matt Ritchie <mritchie@wehi.edu.au>
6
+Author: Benilton S Carvalho <carvalho@bclab.org>, Robert Scharpf <rscharpf@jhsph.edu>, Matt Ritchie <mritchie@wehi.edu.au>, Ingo Ruczinski <iruczins@jhsph.edu>, Rafael A Irizarry
7 7
 Maintainer: Benilton S Carvalho <carvalho@bclab.org>, Robert Scharpf <rscharpf@jhsph.edu>, Matt Ritchie <mritchie@wehi.EDU.AU>
8 8
 Description: Faster implementation of CRLMM specific to SNP 5.0 and 6.0 arrays, as well as a copy number tool specific to 5.0, 6.0, and Illumina platforms
9 9
 License: Artistic-2.0
... ...
@@ -59,11 +59,11 @@ exportClasses(CNSetLM, ffdf, list)
59 59
 exportMethods(open, "[", show, lM, lines, nu, phi, corr, sigma2, tau2)
60 60
 exportMethods(CA, CB, totalCopyNumber)
61 61
 export(crlmm, 
62
-       crlmmCopynumber, 
62
+##       crlmmCopynumber, 
63 63
        crlmmIllumina, 
64 64
        crlmmIllumina2,
65 65
        ellipseCenters,
66
-       genotype, 
66
+##       genotype, 
67 67
        readIdatFiles, 
68 68
        readIdatFiles2,
69 69
        snprma,
... ...
@@ -75,7 +75,7 @@ export(crlmm,
75 75
 export(constructIlluminaCNSet)
76 76
 export(linesCNSetLM)
77 77
 export(computeCN, fit.lm1, fit.lm2, fit.lm3, fit.lm4, construct,
78
-       dqrlsWrapper, nuphiAllele)
78
+       dqrlsWrapper, fit.wls)
79 79
 export(computeCopynumber, ACN)
80 80
 
81 81
 
... ...
@@ -379,11 +379,11 @@ dqrlsWrapper <- function(x, y, wts, tol=1e-7){
379 379
 }
380 380
 
381 381
 fit.wls <- function(allele, Ystar, W, Ns, autosome=TRUE){
382
-	complete <- rowSums(is.na(W)) == 0 
383
-	if(any(!is.finite(W))){## | any(!is.finite(V))){
384
-		i <- which(rowSums(!is.finite(W)) > 0)
385
-		stop("Possible zeros in the within-genotype estimates of the spread (vA, vB). ")
386
-	}
382
+	complete <- which(rowSums(is.na(W)) == 0 & rowSums(is.na(Ystar)) == 0)
383
+##	if(any(!is.finite(W))){## | any(!is.finite(V))){
384
+##		i <- which(rowSums(!is.finite(W)) > 0)
385
+##		stop("Possible zeros in the within-genotype estimates of the spread (vA, vB). ")
386
+##	}
387 387
 	NOHET <- mean(Ns[, 2], na.rm=TRUE) < 0.05
388 388
 	if(missing(allele)) stop("must specify allele")
389 389
 	if(autosome){
... ...
@@ -399,7 +399,7 @@ fit.wls <- function(allele, Ystar, W, Ns, autosome=TRUE){
399 399
 	##How to quickly generate Xstar, Xstar = diag(W) %*% X
400 400
 	##Xstar <- apply(W, 1, generateX, X)
401 401
 	ww <- rep(1, ncol(Ystar))
402
-	for(i in 1:nrow(Ystar)){
402
+	for(i in complete){
403 403
 		betahat[, i] <- dqrlsWrapper(W[i, ] * X, Ystar[i, ], ww)
404 404
 		##ssr <- sum((Ystar[i, ] - matrix(Xstar[, i], nrow(X), ncol(X)) %*% matrix(betahat[, i], ncol(X), 1))^2)
405 405
 	}
... ...
@@ -795,7 +795,7 @@ fit.lm1 <- function(idxBatch,
795 795
 		    MIN.PHI,
796 796
 		    verbose,...){
797 797
 	if(isPackageLoaded("ff")) physical <- get("physical")
798
-	if(verbose) message("Probe batch ", idxBatch, " of ", length(snpBatches))
798
+	if(verbose) message("Probe stratum ", idxBatch, " of ", length(snpBatches))
799 799
 	snps <- snpBatches[[idxBatch]]
800 800
 	batches <- split(seq(along=batch(object)), batch(object))
801 801
 	batches <- batches[sapply(batches, length) >= MIN.SAMPLES]
... ...
@@ -1004,7 +1004,7 @@ fit.lm2 <- function(idxBatch,
1004 1004
 		    MIN.PHI,
1005 1005
 		    verbose,...){
1006 1006
 	physical <- get("physical")
1007
-	if(verbose) message("Probe batch ", idxBatch, " of ", length(snpBatches))
1007
+	if(verbose) message("Probe stratum ", idxBatch, " of ", length(snpBatches))
1008 1008
 	snps <- snpBatches[[idxBatch]]
1009 1009
 	batches <- split(seq(along=batch(object)), batch(object))
1010 1010
 	batches <- batches[sapply(batches, length) >= MIN.SAMPLES]
... ...
@@ -1115,7 +1115,7 @@ fit.lm3 <- function(idxBatch,
1115 1115
 		    MIN.PHI,
1116 1116
 		    verbose, ...){
1117 1117
 	physical <- get("physical")
1118
-	if(verbose) message("Probe batch ", idxBatch, " of ", length(snpBatches))
1118
+	if(verbose) message("Probe stratum ", idxBatch, " of ", length(snpBatches))
1119 1119
 		snps <- snpBatches[[idxBatch]]
1120 1120
 	batches <- split(seq(along=batch(object)), batch(object))
1121 1121
 	batches <- batches[sapply(batches, length) >= MIN.SAMPLES]
... ...
@@ -1354,7 +1354,7 @@ fit.lm4 <- function(idxBatch,
1354 1354
 		    MIN.PHI,
1355 1355
 		    verbose, ...){
1356 1356
 	physical <- get("physical")
1357
-	if(verbose) message("Probe batch ", idxBatch, " of ", length(snpBatches))	
1357
+	if(verbose) message("Probe stratum ", idxBatch, " of ", length(snpBatches))	
1358 1358
 	open(object)
1359 1359
 	open(normal)
1360 1360
 	open(snpflags)
... ...
@@ -1938,48 +1938,48 @@ getFlags <- function(object, PHI.THR){
1938 1938
 ##	return(tmp.objects)
1939 1939
 ##}
1940 1940
 
1941
-##imputeCenter <- function(muA, muB, index.complete, index.missing){
1942
-##	index <- index.missing
1943
-##	mnA <- muA
1944
-##	mnB <- muB
1945
-##	for(j in 1:3){
1946
-##		if(length(index[[j]]) == 0) next()
1947
-##		X <- cbind(1, mnA[index.complete,  -j, drop=FALSE], mnB[index.complete,  -j, drop=FALSE])
1948
-##		Y <- cbind(mnA[index.complete, j], mnB[index.complete,  j])
1949
-##		betahat <- solve(crossprod(X), crossprod(X,Y))
1950
-##		X <- cbind(1, mnA[index[[j]],  -j, drop=FALSE],  mnB[index[[j]],  -j, drop=FALSE])
1951
-##		mus <- X %*% betahat
1952
-##		muA[index[[j]], j] <- mus[, 1]
1953
-##		muB[index[[j]], j] <- mus[, 2]
1954
-##	}
1955
-##	list(muA, muB)
1956
-##}
1957
-##
1958
-##imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
1959
-##	index1 <- which(Ns[, 1] == 0 & Ns[, 3] > MIN.OBS)
1960
-##	if(length(index1) > 0){
1961
-##		X <- cbind(1, muA[index.complete[[1]], 3], muB[index.complete[[1]], 3])
1962
-##		Y <- cbind(1, muA[index.complete[[1]], 1], muB[index.complete[[1]], 1])
1963
-##		betahat <- solve(crossprod(X), crossprod(X,Y))
1964
-##		##now with the incomplete SNPs
1965
-##		X <- cbind(1, muA[index1, 3], muB[index1, 3])
1966
-##		mus <- X %*% betahat
1967
-##		muA[index1, 1] <- mus[, 2]
1968
-##		muB[index1, 1] <- mus[, 3]
1969
-##	}
1970
-##	index1 <- which(Ns[, 3] == 0)
1971
-##	if(length(index1) > 0){
1972
-##		X <- cbind(1, muA[index.complete[[2]], 1], muB[index.complete[[2]], 1])
1973
-##		Y <- cbind(1, muA[index.complete[[2]], 3], muB[index.complete[[2]], 3])
1974
-##		betahat <- solve(crossprod(X), crossprod(X,Y))
1975
-##		##now with the incomplete SNPs
1976
-##		X <- cbind(1, muA[index1, 1], muB[index1, 1])
1977
-##		mus <- X %*% betahat
1978
-##		muA[index1, 3] <- mus[, 2]
1979
-##		muB[index1, 3] <- mus[, 3]
1980
-##	}
1981
-##	list(muA, muB)
1982
-##}
1941
+imputeCenter <- function(muA, muB, index.complete, index.missing){
1942
+	index <- index.missing
1943
+	mnA <- muA
1944
+	mnB <- muB
1945
+	for(j in 1:3){
1946
+		if(length(index[[j]]) == 0) next()
1947
+		X <- cbind(1, mnA[index.complete,  -j, drop=FALSE], mnB[index.complete,  -j, drop=FALSE])
1948
+		Y <- cbind(mnA[index.complete, j], mnB[index.complete,  j])
1949
+		betahat <- solve(crossprod(X), crossprod(X,Y))
1950
+		X <- cbind(1, mnA[index[[j]],  -j, drop=FALSE],  mnB[index[[j]],  -j, drop=FALSE])
1951
+		mus <- X %*% betahat
1952
+		muA[index[[j]], j] <- mus[, 1]
1953
+		muB[index[[j]], j] <- mus[, 2]
1954
+	}
1955
+	list(muA, muB)
1956
+}
1957
+
1958
+imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
1959
+	index1 <- which(Ns[, 1] == 0 & Ns[, 3] > MIN.OBS)
1960
+	if(length(index1) > 0){
1961
+		X <- cbind(1, muA[index.complete[[1]], 3], muB[index.complete[[1]], 3])
1962
+		Y <- cbind(1, muA[index.complete[[1]], 1], muB[index.complete[[1]], 1])
1963
+		betahat <- solve(crossprod(X), crossprod(X,Y))
1964
+		##now with the incomplete SNPs
1965
+		X <- cbind(1, muA[index1, 3], muB[index1, 3])
1966
+		mus <- X %*% betahat
1967
+		muA[index1, 1] <- mus[, 2]
1968
+		muB[index1, 1] <- mus[, 3]
1969
+	}
1970
+	index1 <- which(Ns[, 3] == 0)
1971
+	if(length(index1) > 0){
1972
+		X <- cbind(1, muA[index.complete[[2]], 1], muB[index.complete[[2]], 1])
1973
+		Y <- cbind(1, muA[index.complete[[2]], 3], muB[index.complete[[2]], 3])
1974
+		betahat <- solve(crossprod(X), crossprod(X,Y))
1975
+		##now with the incomplete SNPs
1976
+		X <- cbind(1, muA[index1, 1], muB[index1, 1])
1977
+		mus <- X %*% betahat
1978
+		muA[index1, 3] <- mus[, 2]
1979
+		muB[index1, 3] <- mus[, 3]
1980
+	}
1981
+	list(muA, muB)
1982
+}
1983 1983
 
1984 1984
 ##oneBatch <- function(object, cnOptions, tmp.objects){
1985 1985
 ##	muA <- tmp.objects[["muA"]]
... ...
@@ -119,7 +119,7 @@ not exist. Existence of this file implies that we have already run the
119 119
 copy number estimation and, therefore, we do not need to preprocess
120 120
 and genotype the files a second time.
121 121
 
122
-<<genotype>>=
122
+<<genotype, eval=FALSE>>=
123 123
 if(!file.exists(file.path(outdir, "cnSet.assayData_matrix.rda"))){
124 124
 	gtSet.assayData_matrix <- checkExists("gtSet.assayData_matrix",
125 125
 					   .path=outdir,
... ...
@@ -140,7 +140,7 @@ the files that were on the same chemistry plate.  The code chunk below
140 140
 will load \Robject{cnSet.assayData_matrix} from disk if this
141 141
 computation had already been performed as part of the batch job.
142 142
 
143
-<<copynumber>>=
143
+<<copynumber, eval=FALSE>>=
144 144
 cnSet.assayData_matrix <- checkExists("cnSet.assayData_matrix",
145 145
 				      .path=outdir,
146 146
 				      .FUN=crlmmCopynumber,
... ...
@@ -194,6 +194,7 @@ from the \Robject{outdir}.
194 194
 
195 195
 <<LDS_copynumber>>=
196 196
 Rprof(filename="Rprof_cnff.out", interval=0.1)
197
+##trace(fit.wls, browser)
197 198
 cnSet.assayData_ff <- checkExists("cnSet.assayData_ff",
198 199
 				  .path=outdir,
199 200
 				  .FUN=crlmmCopynumberLD,