Browse code

Exporting crlmmIllumina2. Added constructIlluminaCNSet function.

illumina_copynumber vignette is using the functions readIdatFiles2 and crlmmIllumina2.

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

Rob Scharp authored on 02/08/2010 08:52:15
Showing3 changed files

... ...
@@ -60,6 +60,7 @@ exportMethods(open, "[", show, lM, lines)
60 60
 export(crlmm, 
61 61
        crlmmCopynumber, 
62 62
        crlmmIllumina, 
63
+       crlmmIllumina2,
63 64
        genotype, 
64 65
        readIdatFiles, 
65 66
        readIdatFiles2,
... ...
@@ -69,6 +70,7 @@ export(crlmm,
69 70
        genotype2, genotypeLD,
70 71
        batch,
71 72
        crlmmCopynumber2, crlmmCopynumberLD)
73
+export(constructIlluminaCNSet)
72 74
 
73 75
 
74 76
 
... ...
@@ -3159,3 +3159,23 @@ constructIlluminaAssayData <- function(np, snp, object, storage.mode="environmen
3159 3159
 			   CA=emptyMatrix,
3160 3160
 			   CB=emptyMatrix)
3161 3161
 }
3162
+constructIlluminaCNSet <- function(res,
3163
+				   cnAB,
3164
+				   crlmmResult){
3165
+	fD <- constructIlluminaFeatureData(c(res$gns, cnAB$gns), cdfName="human370v1c")
3166
+	new.order <- order(fD$chromosome, fD$position)
3167
+	fD <- fD[new.order, ]
3168
+	aD <- constructIlluminaAssayData(cnAB, res, crlmmResult, order.index=new.order)
3169
+	protocolData(crlmmResult)$batch <- vector("integer", ncol(crlmmResult))
3170
+	container <- new("CNSetLM", 
3171
+			 assayData=aD,
3172
+			 phenoData=phenoData(crlmmResult),
3173
+			 protocolData=protocolData(crlmmResult),
3174
+			 featureData=fD,
3175
+			 annotation="human370v1c")
3176
+	lM(container) <- initializeParamObject(list(featureNames(container), unique(protocolData(container)$batch)))
3177
+	container
3178
+}
3179
+				   
3180
+				   
3181
+
... ...
@@ -46,6 +46,7 @@ options(continue=" ")
46 46
 @ 
47 47
 
48 48
 <<libraries>>=
49
+library(ff)
49 50
 library(crlmm)
50 51
 crlmm:::validCdfNames()
51 52
 if(length(grep("development", sessionInfo()[[1]]$status)) == 1){
... ...
@@ -92,7 +93,7 @@ redfiles = all(file.exists(paste(arrayNames, "_Red.idat", sep="")))
92 93
 <<samplesToProcess2, echo=FALSE>>=
93 94
 ## To speed up repeated calls to Sweave 
94 95
 RG <- checkExists("RG", .path=outdir,
95
-		  .FUN=readIdatFiles,
96
+		  .FUN=readIdatFiles2,
96 97
 		  sampleSheet=samplesheet,
97 98
 		  path=dirname(arrayNames[1]),
98 99
 		  arrayInfoColNames=list(barcode=NULL, position="SentrixPosition"), 
... ...
@@ -100,7 +101,7 @@ RG <- checkExists("RG", .path=outdir,
100 101
 annotation(RG) <- "human370v1c"
101 102
 crlmmResult <- checkExists("crlmmResult",
102 103
 			   .path=outdir,
103
-			   .FUN=crlmmIllumina,
104
+			   .FUN=crlmmIllumina2,
104 105
 			   RG=RG,
105 106
 			   sns=pData(RG)$ID, 
106 107
 			   returnParams=TRUE,
... ...
@@ -111,7 +112,7 @@ protocolData(crlmmResult)$ScanDate <- protocolData(RG)$ScanDate
111 112
 @ 
112 113
 
113 114
 <<samplesToProcess3, eval=FALSE>>=
114
-RG <- readIdatFiles(samplesheet, 
115
+RG <- readIdatFiles2(samplesheet[1:10, ], 
115 116
 		    path=dirname(arrayNames[1]), 
116 117
 		    arrayInfoColNames=list(barcode=NULL, position="SentrixPosition"), 
117 118
 		    saveDate=TRUE)
... ...
@@ -136,8 +137,9 @@ cnAB <- get("cnAB")
136 137
 @ 
137 138
 
138 139
 <<loadIntermediate, echo=FALSE>>=
139
-res <- checkExists("res", .path=outdir, .FUN=load, file=file.path(outdir, "snpFile.rda"))
140
-cnAB <- checkExists("cnAB", .path=outdir, .FUN=load, file=file.path(outdir, "cnFile.rda"))
140
+trace(checkExists, browser)
141
+res <- checkExists("snpFile", .path=outdir, .FUN=load, file=file.path(outdir, "snpFile.rda"))
142
+cnAB <- checkExists("cnFile", .path=outdir, .FUN=load, file=file.path(outdir, "cnFile.rda"))
141 143
 @ 
142 144
 
143 145
 After running the crlmm algorithm, we construct a container for
... ...
@@ -148,42 +150,19 @@ inst/scripts directory of this package and can be sourced as follows.
148 150
 Documentation for these helper functions will be available in the
149 151
 devel version of this package.
150 152
 
151
-<<constructContainer, eval=FALSE>>=
152
-path <- system.file("scripts", package="crlmm")
153
-source(file.path(path, "helperFunctions.R"))
154
-fD <- constructFeatureData(c(res$gns, cnAB$gns), cdfName="human370v1c")
155
-new.order <- order(fD$chromosome, fD$position)
156
-fD <- fD[new.order, ]
157
-aD <- constructAssayData(cnAB, res, crlmmResult, order.index=new.order)
158
-protocolData(crlmmResult)$batch <- vector("integer", ncol(crlmmResult))
159
-container <- new("CNSetLM", 
160
-		 assayData=aD,
161
-		 phenoData=phenoData(crlmmResult),
162
-		 protocolData=protocolData(crlmmResult),
163
-		 featureData=fD,
164
-		 annotation="human370v1c")
153
+<<constructContainer1, check=FALSE>>=
154
+container <- checkExists("container",
155
+			 .path=outdir,
156
+			 .FUN=constructIlluminaCNSet,
157
+			 res=res,
158
+			 cnAB=cnAB,
159
+			 crlmmResult=crlmmResult)
165 160
 @ 
166 161
 
167
-<<container2,echo=FALSE>>=
168
-if(!file.exists(file.path(outdir, "cnSet"))){
169
-	path <- system.file("scripts", package="crlmm")
170
-	source(file.path(path, "helperFunctions.R"))
171
-	fD <- constructFeatureData(c(res$gns, cnAB$gns), cdfName="human370v1c")
172
-	new.order <- order(fD$chromosome, fD$position)
173
-	fD <- fD[new.order, ]
174
-	aD <- constructAssayData(cnAB, res, crlmmResult, order.index=new.order)
175
-	protocolData(crlmmResult)$batch <- vector("integer", ncol(crlmmResult))
176
-	container <- new("CNSetLM", 
177
-			 assayData=aD,
178
-			 phenoData=phenoData(crlmmResult),
179
-			 protocolData=protocolData(crlmmResult),
180
-			 featureData=fD,
181
-			 annotation="human370v1c")
182
-}
162
+<<constructContainer2, eval=FALSE>>=
163
+container <- constructIlluminaCNSet(res=res, cnAB=cnAB, crlmmResult=crlmmResult)
183 164
 @ 
184 165
 
185
-
186
-
187 166
 As described in the \texttt{copynumber} vignette, two \R{} functions
188 167
 for copy number estimation are available: \Rfunction{crlmmCopynumber}
189 168
 and \Rfunction{crlmmCopynumber2}.  The latter requires that the assay
... ...
@@ -196,17 +175,14 @@ decided at the beginning of the analysis and then propogated to both
196 175
 the genotyping and copy number estimation steps.  Here, we use the
197 176
 \Rfunction{crlmmCopynumber} to estimate copy number.
198 177
 
199
-<<estimateCopynumber, echo=FALSE>>=
200
-if(!exists("cnSet")){
201
-	if(!file.exists(file.path(outdir, "cnSet.rda"))){
202
-		cnSet <- crlmmCopynumber(container, verbose=FALSE)
203
-		save(cnSet, file=file.path(outdir, "cnSet.rda"))
204
-	} else load(file.path(outdir, "cnSet.rda"))
205
-}
178
+<<estimateCopynumber1, echo=FALSE>>=
179
+cnSet <- checkExists("cnSet", .path=outdir,
180
+		     .FUN=crlmmCopynumberLD,
181
+		     object=container)
206 182
 @ 
207 183
 
208 184
 <<estimateCopynumber, eval=FALSE>>=
209
-cnSet <- crlmmCopynumber(container, verbose=TRUE)
185
+cnSet <- crlmmCopynumberLD(container)
210 186
 @ 
211 187
 
212 188
 %In the following code, we define helper functions to construct a