Browse code

removed weighted.lm option for crlmmCopynumber

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

Rob Scharp authored on 21/08/2010 02:48:39
Showing 2 changed files

... ...
@@ -657,8 +657,7 @@ crlmmCopynumberLD <- function(object,
657 657
 			    MIN.NU=2^3,
658 658
 			    MIN.PHI=2^3,
659 659
 			    THR.NU.PHI=TRUE,
660
-			    thresholdCopynumber=TRUE,
661
-			      weighted.lm=TRUE){
660
+			    thresholdCopynumber=TRUE){
662 661
 	stopifnot("batch" %in% varLabels(protocolData(object)))
663 662
 	stopifnot("chromosome" %in% fvarLabels(object))
664 663
 	stopifnot("position" %in% fvarLabels(object))
... ...
@@ -794,8 +793,7 @@ fit.lm1 <- function(idxBatch,
794 793
 		    THR.NU.PHI,
795 794
 		    MIN.NU,
796 795
 		    MIN.PHI,
797
-		    verbose,
798
-		    weighted.lm, ...){
796
+		    verbose,...){
799 797
 	if(isPackageLoaded("ff")) physical <- get("physical")
800 798
 	if(verbose) message("Probe batch ", idxBatch, " of ", length(snpBatches))
801 799
 	snps <- snpBatches[[idxBatch]]
... ...
@@ -907,20 +905,18 @@ fit.lm1 <- function(idxBatch,
907 905
 		wB <- sqrt(1/vB2)
908 906
 		YA <- muA*wA
909 907
 		YB <- muB*wB
910
-		if(weighted.lm){
911
-			res <- fit.wls(allele="A", Ystar=YA, W=wA, Ns=Ns)
912
-		} else{
913
-			if(zzzz==1) message("currently, only weighted least squares (wls) is available... fitting wls")
914
-			res <- fit.wls(allele="A", Ystar=YA, W=wA, Ns=Ns)			
915
-		}
908
+		res <- fit.wls(allele="A", Ystar=YA, W=wA, Ns=Ns)
909
+##		} else{
910
+##			if(zzzz==1) message("currently, only weighted least squares (wls) is available... fitting wls")
911
+##			res <- fit.wls(allele="A", Ystar=YA, W=wA, Ns=Ns)			
912
+##		}
916 913
 		nuA[, J] <- res[[1]]
917 914
 		phiA[, J] <- res[[2]]
918
-		if(weighted.lm){
919
-			res <- fit.wls(allele="B", Ystar=YB, W=wB, Ns=Ns)
920
-		} else {
921
-			if(zzzz==1) message("currently, only weighted least squares (wls) is available... fitting wls")
922
-			res <- fit.wls(allele="B", Ystar=YB, W=wB, Ns=Ns)			
923
-		}
915
+		res <- fit.wls(allele="B", Ystar=YB, W=wB, Ns=Ns)
916
+##		} else {
917
+##			if(zzzz==1) message("currently, only weighted least squares (wls) is available... fitting wls")
918
+##			res <- fit.wls(allele="B", Ystar=YB, W=wB, Ns=Ns)			
919
+##		}
924 920
 		##nuB[, J] <- res[[1]]
925 921
 		nuB[, J] <- res[1, ]
926 922
 		##phiB[, J] <- res[[2]]
... ...
@@ -1006,8 +1002,7 @@ fit.lm2 <- function(idxBatch,
1006 1002
 		    THR.NU.PHI,
1007 1003
 		    MIN.NU,
1008 1004
 		    MIN.PHI,
1009
-		    verbose,
1010
-		    weighted.lm, ...){
1005
+		    verbose,...){
1011 1006
 	physical <- get("physical")
1012 1007
 	if(verbose) message("Probe batch ", idxBatch, " of ", length(snpBatches))
1013 1008
 	snps <- snpBatches[[idxBatch]]
... ...
@@ -2815,8 +2810,7 @@ computeCN <- function(object,
2815 2810
 		      MIN.PHI=2^3,
2816 2811
 		      THR.NU.PHI=TRUE,
2817 2812
 		      thresholdCopynumber=TRUE,
2818
-		      type=c("autosome.snps", "autosome.nps", "X.snps", "X.nps"),
2819
-		      weighted.lm=TRUE){
2813
+		      type=c("autosome.snps", "autosome.nps", "X.snps", "X.nps")){
2820 2814
 	stopifnot("batch" %in% varLabels(protocolData(object)))
2821 2815
 	stopifnot("chromosome" %in% fvarLabels(object))
2822 2816
 	stopifnot("position" %in% fvarLabels(object))
... ...
@@ -2930,8 +2924,7 @@ computeCN <- function(object,
2930 2924
 		 MIN.NU=MIN.NU,
2931 2925
 		 MIN.PHI=MIN.PHI,
2932 2926
 		 verbose=verbose,
2933
-		 neededPkgs="crlmm",
2934
-		 weighted.lm=weighted.lm)
2927
+		 neededPkgs="crlmm")
2935 2928
 	message("finished")
2936 2929
 	return(obj)
2937 2930
 }
... ...
@@ -120,7 +120,7 @@ copy number estimation and, therefore, we do not need to preprocess
120 120
 and genotype the files a second time.
121 121
 
122 122
 <<genotype>>=
123
-if(!file.exists(file.path(outdir, "cnSet.rda"))){
123
+if(!file.exists(file.path(outdir, "cnSet.assayData_matrix.rda"))){
124 124
 	gtSet.assayData_matrix <- checkExists("gtSet.assayData_matrix",
125 125
 					   .path=outdir,
126 126
 					   .FUN=genotype,
... ...
@@ -146,10 +146,6 @@ cnSet.assayData_matrix <- checkExists("cnSet.assayData_matrix",
146 146
 				      .FUN=crlmmCopynumber,
147 147
 				      object=gtSet.assayData_matrix,
148 148
 				      chromosome=22)
149
-##Rprof(interval=0.1)
150
-obj <- crlmmCopynumber(gtSet.assayData_matrix, chromosome=23)
151
-if(file.exists(file.path(outdir, "gtSet.assayData_matrix.rda")))
152
-	unlink(file.path(outdir, "gtSet.assayData_matrix.rda"))
153 149
 @ 
154 150
 
155 151
 
... ...
@@ -201,13 +197,12 @@ Rprof(filename="Rprof_cnff.out", interval=0.1)
201 197
 cnSet.assayData_ff <- checkExists("cnSet.assayData_ff",
202 198
 				  .path=outdir,
203 199
 				  .FUN=crlmmCopynumberLD,
204
-				  filenames=celFiles,
205
-				  cdfName=cdfName,
206
-				  batch=batch)
200
+				  object=gtSet.assayData_ff)
207 201
 Rprof(NULL)
208
-if(file.exists(file.path(outdir, "gtSet.assayData_ff.rda")))
209
-	unlink(file.path(outdir, "gtSet.assayData_ff.rda"))
202
+##if(file.exists(file.path(outdir, "gtSet.assayData_ff.rda")))
203
+##	unlink(file.path(outdir, "gtSet.assayData_ff.rda"))
210 204
 gc()
205
+q("no")
211 206
 @ 
212 207
 
213 208
 The objects returned by the \Rfunction{genotypeLD} and