Browse code

commented CA<- CB<- in nonpolymorphic, polymorphic, and computeCopynumber functions

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

Rob Scharp authored on 21/08/2010 02:48:04
Showing 3 changed files

... ...
@@ -57,7 +57,7 @@ construct <- function(filenames, cdfName, copynumber=FALSE,
57 57
 		stopifnot(length(batch) == length(sns))
58 58
 	}
59 59
 	if(missing(sns) & missing(filenames)) stop("one of filenames or samplenames (sns) must be provided")
60
-	if(verbose) message("Initializing container for assay data elements alleleA, alleleB, call, callProbability, CA, CB")
60
+	if(verbose) message("Initializing container for assay data elements alleleA, alleleB, call, callProbability")
61 61
 	if(!missing(filenames)){
62 62
 		if(missing(sns)) sns <- basename(filenames)
63 63
 		protocolData <- getProtocolData.Affy(filenames)
... ...
@@ -605,15 +605,15 @@ crlmmCopynumber <- function(object,
605 605
 			if(verbose) message("Batch ", ii, " of ", length(which.batches))
606 606
 			row.index <- which(chromosome(object) == i)
607 607
 			##Note that ffdf assayDataElements are data.frames after subsetting(not matrices)
608
-			ca <- as.matrix(CA(object)[row.index, j])
609
-			cb <- as.matrix(CB(object)[row.index, j])
610
-			dimnames(ca) <- dimnames(cb) <- list(featureNames(object)[row.index], sampleNames(object)[j])
608
+##			ca <- as.matrix(CA(object)[row.index, j])
609
+##			cb <- as.matrix(CB(object)[row.index, j])
610
+##			dimnames(ca) <- dimnames(cb) <- list(featureNames(object)[row.index], sampleNames(object)[j])
611 611
 			tmp <- new("CNSet",
612 612
 				   call=as.matrix(calls(object)[row.index, j]),
613 613
 				   callProbability=as.matrix(snpCallProbability(object)[row.index, j]),
614 614
 				   alleleA=as.matrix(A(object)[row.index, j]),
615 615
 				   alleleB=as.matrix(B(object)[row.index, j]),
616
-				   CA=ca, CB=cb,
616
+##				   CA=ca, CB=cb,
617 617
 				   phenoData=phenoData(object)[j, ],
618 618
 				   annotation=annotation(object))
619 619
 			featureData(tmp) <- addFeatureAnnotation(tmp)
... ...
@@ -634,12 +634,12 @@ crlmmCopynumber <- function(object,
634 634
 						 THR.NU.PHI=THR.NU.PHI,
635 635
 						 thresholdCopynumber=thresholdCopynumber)
636 636
 			fData(tmp) <- fData(tmp)[, -(1:3)]
637
-			CA(tmp) <- matrix(as.integer(CA(tmp)*100), nrow=nrow(tmp), ncol=ncol(tmp),
638
-					  dimnames=list(featureNames(tmp), sampleNames(tmp)))
639
-			CB(tmp) <- matrix(as.integer(CB(tmp)*100), nrow=nrow(tmp), ncol=ncol(tmp),
640
-					  dimnames=list(featureNames(tmp), sampleNames(tmp)))
641
-			CA(object)[row.index, j] <- CA(tmp)
642
-			CB(object)[row.index, j] <- CB(tmp)
637
+##			CA(tmp) <- matrix(as.integer(CA(tmp)*100), nrow=nrow(tmp), ncol=ncol(tmp),
638
+##					  dimnames=list(featureNames(tmp), sampleNames(tmp)))
639
+##			CB(tmp) <- matrix(as.integer(CB(tmp)*100), nrow=nrow(tmp), ncol=ncol(tmp),
640
+##					  dimnames=list(featureNames(tmp), sampleNames(tmp)))
641
+##			CA(object)[row.index, j] <- CA(tmp)
642
+##			CB(object)[row.index, j] <- CB(tmp)
643 643
 			##ad-hocery
644 644
 			batchName <- unique(batch(object)[j])
645 645
 			fvarLabels(tmp)[15:17] <- paste(c("corrAB", "corrBB", "corrAA"), batchName, sep=".")
... ...
@@ -1744,8 +1744,8 @@ nonpolymorphic <- function(object, cnOptions, tmp.objects){
1744 1744
 	muB <- tmp.objects[["muB"]]
1745 1745
 	A <- A(object)
1746 1746
 	B <- B(object)
1747
-	CA <- CA(object)
1748
-	CB <- CB(object)
1747
+##	CA <- CA(object)
1748
+##	CB <- CB(object)
1749 1749
 	if(CHR == 23){
1750 1750
 		phiAX <- getParam(object, "phiAX", batch)
1751 1751
 		phiBX <- getParam(object, "phiBX", batch)
... ...
@@ -1813,20 +1813,12 @@ nonpolymorphic <- function(object, cnOptions, tmp.objects){
1813 1813
 		if(any(pseudoAR)){
1814 1814
 			nu1[pseudoAR] <- 2^(mus[pseudoAR, 1]) - 2*phi1[pseudoAR]
1815 1815
 		}
1816
-		CT1 <- 1/phi1*(A.male-nu1)
1817
-		CT2 <- 1/phi2*(A.female-nu2)
1818
-		##CT2 <- 1/phi2*(NP[, gender=="female"]-nu2)
1819
-		##CT1 <- matrix(as.integer(100*CT1), nrow(CT1), ncol(CT1))
1820
-		##CT2 <- matrix(as.integer(100*CT2), nrow(CT2), ncol(CT2))
1821
-		##CT <- envir[["CT"]]
1822
-		CA <- CA(obj1)
1823
-		CA[, gender==1] <- CT1
1824
-		CA[, gender==2] <- CT2
1825
-		CA(object)[!isSnp(object), ] <- CA
1826
-		##CT[, plate==uplate[p] & gender=="male"] <- CT1
1827
-		##CT[, plate==uplate[p] & gender=="female"] <- CT2
1828
-		##envir[["CT"]] <- CT
1829
-
1816
+##		CT1 <- 1/phi1*(A.male-nu1)
1817
+##		CT2 <- 1/phi2*(A.female-nu2)
1818
+##		CA <- CA(obj1)
1819
+##		CA[, gender==1] <- CT1
1820
+##		CA[, gender==2] <- CT2
1821
+##		CA(object)[!isSnp(object), ] <- CA
1830 1822
 		##only using females to compute the variance
1831 1823
 		##normalNP[, gender=="male"] <- NA
1832 1824
 		normal[, gender==1] <- NA
... ...
@@ -2287,13 +2279,13 @@ polymorphic <- function(object, cnOptions, tmp.objects){
2287 2279
 	
2288 2280
 	nuA <- getParam(object, "nuA", batch)
2289 2281
 	nuB <- getParam(object, "nuB", batch)
2290
-	nuA.se <- getParam(object, "nuA.se", batch)
2291
-	nuB.se <- getParam(object, "nuB.se", batch)
2282
+##	nuA.se <- getParam(object, "nuA.se", batch)
2283
+##	nuB.se <- getParam(object, "nuB.se", batch)
2292 2284
 
2293 2285
 	phiA <- getParam(object, "phiA", batch)
2294 2286
 	phiB <- getParam(object, "phiB", batch)
2295
-	phiA.se <- getParam(object, "phiA.se", batch)
2296
-	phiB.se <- getParam(object, "phiB.se", batch)
2287
+##	phiA.se <- getParam(object, "phiA.se", batch)
2288
+##	phiB.se <- getParam(object, "phiB.se", batch)
2297 2289
 	A <- A(object)
2298 2290
 	B <- B(object)
2299 2291
 
... ...
@@ -2305,14 +2297,14 @@ polymorphic <- function(object, cnOptions, tmp.objects){
2305 2297
 		phiAX <- getParam(object, "phiAX", batch)  ##nonspecific hybridization coef
2306 2298
 		phiBX <- getParam(object, "phiBX", batch)  ##nonspecific hybridization coef			
2307 2299
 		phistar <- phiBX/phiA  
2308
-		tmp <- (B-nuB - phistar*A + phistar*nuA)/phiB
2309
-		copyB <- tmp/(1-phistar*phiAX/phiB)
2310
-		copyA <- (A-nuA-phiAX*copyB)/phiA
2311
-		CB(object) <- copyB  ## multiplies by 100 and converts to integer 
2312
-		CA(object) <- copyA
2300
+##		tmp <- (B-nuB - phistar*A + phistar*nuA)/phiB
2301
+##		copyB <- tmp/(1-phistar*phiAX/phiB)
2302
+##		copyA <- (A-nuA-phiAX*copyB)/phiA
2303
+##		CB(object) <- copyB  ## multiplies by 100 and converts to integer 
2304
+##		CA(object) <- copyA
2313 2305
 	} else{
2314
-		CA(object) <- matrix((1/phiA*(A-nuA)), nrow(A), ncol(A))
2315
-		CB(object) <- matrix((1/phiB*(B-nuB)), nrow(B), ncol(B))
2306
+##		CA(object) <- matrix((1/phiA*(A-nuA)), nrow(A), ncol(A))
2307
+##		CB(object) <- matrix((1/phiB*(B-nuB)), nrow(B), ncol(B))
2316 2308
 
2317 2309
 	}
2318 2310
 	return(object)
... ...
@@ -2667,28 +2659,7 @@ thresholdModelParams <- function(object, cnOptions){
2667 2659
 	return(object)
2668 2660
 }
2669 2661
 
2670
-##computeCopynumber.SnpSuperSet <- function(object, cnOptions){
2671
-####	use.ff <- cnOptions[["use.ff"]]
2672
-####	if(!use.ff){
2673
-####		object <- as(object, "CrlmmSet")
2674
-####	} else	object <- as(object, "CrlmmSetFF")
2675
-##	bias.adj <- cnOptions[["bias.adj"]]
2676
-##	##must be FALSE to initialize parameters
2677
-##	cnOptions[["bias.adj"]] <- FALSE
2678
-##	## Add linear model parameters to the CrlmmSet object
2679
-##	featureData(object) <- lm.parameters(object, cnOptions)
2680
-##	if(!isValidCdfName(annotation(object))) stop(annotation(object), " not supported.")
2681
-##	object <- as(object, "CNSet")
2682
-##	object <- computeCopynumber.CNSet(object, cnOptions)
2683
-##	if(bias.adj==TRUE){## run a second time
2684
-##		object <- computeCopynumber.CNSet(object, cnOptions)
2685
-##	}
2686
-##	return(object)
2687
-##}
2688
-
2689
-
2690 2662
 computeCopynumber.CNSet <- function(object, cnOptions){
2691
-	##PLATE <- unique(object$batch)
2692 2663
 	PLATE <- unique(batch(object))
2693 2664
 	verbose <- cnOptions$verbose
2694 2665
 	tmp.objects <- instantiateObjects(object, cnOptions)
... ...
@@ -93,7 +93,7 @@ setMethod("computeCopynumber", "CNSet",
93 93
 				    THR.NU.PHI=THR.NU.PHI,
94 94
 				    thresholdCopynumber=thresholdCopynumber)
95 95
 	bias.adj <- cnOptions[["bias.adj"]]
96
-	if(bias.adj & all(is.na(CA(object)))){
96
+	if(bias.adj & all(is.na(nu(object, "A")[, 1])){
97 97
 		cnOptions[["bias.adj"]] <- FALSE
98 98
 	}
99 99
 	object <- computeCopynumber.CNSet(object, cnOptions)				
... ...
@@ -129,7 +129,8 @@ if(!file.exists(file.path(outdir, "cnSet.rda"))){
129 129
 					   batch=batch[1:20])
130 130
 	class(calls(gtSet.assayData_matrix))
131 131
 }
132
-q("no")
132
+##q("no")
133
+stop()
133 134
 @ 
134 135
 
135 136
 Next, we estimate copy number for the 20 CEL files.  The copy number