Browse code

CB estimates were NA for SNP indices were mixed with NP indices. Bug was fixed in the ACN function

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

Rob Scharp authored on 30/08/2010 19:40:54
Showing6 changed files

... ...
@@ -87,7 +87,12 @@ exportMethods(A, B, nuA, nuB, phiA, phiB, corr, tau2, Ns, medians, mads)
87 87
 export(genotypeSummary,
88 88
        shrinkSummary,
89 89
        estimateCnParameters,
90
-       shrinkGenotypeSummaries, indexComplete, summarizeSnps, constructIlluminaAssayData, crlmmIlluminaV2, ACN, C1, C2, C3, fit.lm3)
90
+       shrinkGenotypeSummaries,
91
+##       indexComplete,
92
+       summarizeSnps,
93
+       constructIlluminaAssayData,
94
+       ACN, C1, C2, C3, fit.lm3,
95
+       linesCNSet)
91 96
 
92 97
 exportMethods(N.AA, N.AB, N.BB, "N.AA<-", "N.AB<-", "N.BB<-", medianA.AA, "medianA.AA<-")
93 98
 
... ...
@@ -43,7 +43,7 @@ setMethod("Ns", signature(object="AssayData"),
43 43
 		  return(res)
44 44
 	  })
45 45
 setMethod("corr", signature(object="AssayData"),
46
-	  function(object, i, j, ...){
46
+	  function(object, ...){
47 47
 		  dotArgs <- names(list(...))
48 48
 		  missing.i <- !("i" %in% dotArgs)
49 49
 		  missing.j <- !("j" %in% dotArgs)
... ...
@@ -410,7 +410,6 @@ ACN <- function(object, allele, i , j){
410 410
 		is.X <- chromosome(object)[i]==23 & is.ann
411 411
 		is.auto <- chromosome(object)[i] < 23 & is.ann
412 412
 		is.snp <- isSnp(object)[i] & is.ann
413
-
414 413
 	} else{
415 414
 		is.ann <- !is.na(chromosome(object))
416 415
 		is.X <- chromosome(object)==23 & is.ann
... ...
@@ -509,8 +508,8 @@ ACN <- function(object, allele, i , j){
509 508
 			acn.index <- which(!is.snp)
510 509
 			marker.index <- i[!is.snp]
511 510
 			acn[acn.index, ] <- 0
512
-		} else{
513
-			## SNPs
511
+		}
512
+		if(any(is.snp)){
514 513
 			if(any(is.auto)){
515 514
 				## autosomal SNPs
516 515
 				acn.index <- which(is.auto & is.snp)
... ...
@@ -6,21 +6,28 @@ setMethod("lines", signature=signature(x="CNSet"),
6 6
 linesCNSet <- function(x, y, batch, copynumber, x.axis="A", ...){
7 7
 	require(ellipse)
8 8
 	object <- x
9
-	I <- y
10
-	stopifnot(length(I) == 1)
11
-	column <- grep(batch, batchNames(object))
12
-	stopifnot(length(column) == 1)
13
-	nuA <- nu(object, "A")[I, column]
14
-	nuB <- nu(object, "B")[I, column]
15
-	phiA <- phi(object, "A")[I, column]
16
-	phiB <- phi(object, "B")[I, column]
17
-	tau2A <- tau2(object, "A")[I, column]
18
-	tau2B <- tau2(object, "B")[I, column]
19
-	sigma2A <- sigma2(object, "A")[I, column]
20
-	sigma2B <- sigma2(object, "B")[I, column]
21
-	corrAB <- corr(object, "AB")[I, column]
22
-	corrAA <- corr(object, "AA")[I, column]
23
-	corrBB <- corr(object, "BB")[I, column]
9
+	marker.index <- y
10
+	stopifnot(length(marker.index) == 1)
11
+	batch.index <- match(batch, batchNames(object))
12
+	stopifnot(length(batch.index) == 1)
13
+	nuA <- nu(object, "A")[marker.index, batch.index]
14
+	nuB <- nu(object, "B")[marker.index, batch.index]
15
+	phiA <- phi(object, "A")[marker.index, batch.index]
16
+	phiB <- phi(object, "B")[marker.index, batch.index]
17
+
18
+	taus <- tau2(object, i=marker.index, j=batch.index)
19
+	tau2A <- taus[, "A", "BB", ]
20
+	tau2B <- taus[, "B", "AA", ]
21
+	sigma2A <- taus[, "A", "AA", ]
22
+	sigma2B <- taus[, "B", "BB", ]
23
+##	tau2A <- tau2(object, "A")[marker.index, batch.index]
24
+##	tau2B <- tau2(object, "B")[marker.index, batch.index]
25
+##	sigma2A <- sigma2(object, "A")[marker.index, batch.index]
26
+##	sigma2B <- sigma2(object, "B")[marker.index, batch.index]
27
+	cors <- corr(object, i=marker.index, j=batch.index)[, , ]
28
+	corrAB <- cors[["AB"]]
29
+	corrAA <- cors[["AA"]]
30
+	corrBB <- cors[["BB"]]
24 31
 	if(all(is.na(nuA))) {
25 32
 		message("Parameter estimates for batch ", batch, " not available")
26 33
 		next()
... ...
@@ -47,7 +47,7 @@ options(continue=" ")
47 47
 @
48 48
 
49 49
 <<libraries>>=
50
-library(ff)
50
+##library(ff)
51 51
 library(crlmm)
52 52
 crlmm:::validCdfNames()
53 53
 if(getRversion() < "2.12.0"){
... ...
@@ -57,7 +57,7 @@ outdir <- paste("/thumper/ctsa/snpmicroarray/rs/ProcessedData/crlmm/", rpath, "/
57 57
 datadir <- "/thumper/ctsa/snpmicroarray/illumina/IDATS/370k"
58 58
 @
59 59
 
60
-Options for controlling RAM with the \Rpackage{ff} package.
60
+%Options for controlling RAM with the \Rpackage{ff} package.
61 61
 
62 62
 <<ldOptions>>=
63 63
 ldPath(outdir)
... ...
@@ -114,18 +114,41 @@ grnfiles = all(file.exists(paste(arrayNames, "_Grn.idat", sep="")))
114 114
 redfiles = all(file.exists(paste(arrayNames, "_Red.idat", sep="")))
115 115
 @
116 116
 
117
-<<samplesToProcess2>>=
117
+<<eval=FALSE>>=
118
+RG <- checkExists("RG", .path=outdir,
119
+		  .FUN=readIdatFiles,
120
+		  sampleSheet=samplesheet,
121
+		  path=dirname(arrayNames[1]),
122
+		  arrayInfoColNames=list(barcode=NULL, position="SentrixPosition"),
123
+		  saveDate=TRUE)
124
+annotation(RG) <- "human370v1c"
118 125
 crlmmResult <- checkExists("crlmmResult",
119 126
 			   .path=outdir,
120
-			   .FUN=crlmm:::crlmmIlluminaV2,
127
+			   .FUN=crlmmIllumina,
128
+			   RG=RG,
129
+			   sns=pData(RG)$ID,
130
+			   returnParams=TRUE,
131
+			   cnFile=file.path(outdir, "cnFile.rda"),
132
+			   snpFile=file.path(outdir, "snpFile.rda"),
133
+			   save.it=TRUE)
134
+protocolData(crlmmResult)$ScanDate <- protocolData(RG)$ScanDate
135
+@
136
+
137
+<<samplesToProcess2, eval=TRUE, echo=FALSE>>=
138
+crlmmResult <- checkExists("crlmmResult",
139
+			   .path=outdir,
140
+			   .FUN=crlmmIlluminaV2,
121 141
 			   sampleSheet=samplesheet,
122 142
 			   path=dirname(arrayNames[1]),
123 143
 			   arrayInfoColNames=list(barcode=NULL, position="SentrixPosition"),
124 144
 			   returnParams=TRUE,
125 145
 			   cnFile=file.path(outdir, "cnFile.rda"),
126 146
 			   snpFile=file.path(outdir, "snpFile.rda"),
127
-			   save.ab=TRUE,
147
+			   save.it=TRUE,
128 148
 			   cdfName="human370v1c", .load.it=FALSE)
149
+## MR: shouldn't this be intensities in the A element of res?
150
+load(file.path(outdir, "snpFile.rda"))
151
+res[["A"]][1:5, 1:5]
129 152
 stop()
130 153
 @
131 154
 
... ...
@@ -151,9 +174,10 @@ facilitating the construction of this container have been added to the
151 174
 inst/scripts directory of this package and can be sourced as follows.
152 175
 
153 176
 <<constructContainer2>>=
177
+load(file.path(outdir, "snpFile.rda"))
154 178
 batch <- as.factor(rep("1", ncol(crlmmResult)))
155 179
 container <- constructIlluminaCNSet(crlmmResult=crlmmResult, snpFile=file.path(outdir, "snpFile.rda"), cnFile=file.path(outdir, "cnFile.rda"), batch=batch)
156
-stop()
180
+A(container)[1:5, 1:5]
157 181
 @
158 182
 
159 183
 As described in the \texttt{copynumber} vignette, two \R{} functions
... ...
@@ -30,7 +30,9 @@ data(sample.CNSetLM)
30 30
 ## update to class CNSet
31 31
 cnSet <- as(sample.CNSetLM, "CNSet")
32 32
 all(isCurrent(cnSet)) ## is the cnSet object current?
33
-
33
+##subsetting
34
+cnSet2 <- cnSet[, 1:5]
35
+stopifnot(batchNames(cnSet2) == "C")
34 36
 \dontrun{
35 37
 	## updating class CNSetLM using ff objects
36 38
 	## a bigger object with multiple batches