Browse code

Added constructIlluminaFeatureData and constructIlluminaAssayData to cnrma-functions

A few edits to illumina_copynumber vignette:
- checkExists for RG
- checkExists for crlmmResult
- checkExists for loading snpFile and cnFile

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

Rob Scharp authored on 02/08/2010 08:52:11
Showing 4 changed files

... ...
@@ -1,7 +1,7 @@
1 1
 Package: crlmm
2 2
 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
-Version: 1.7.8
4
+Version: 1.7.9
5 5
 Date: 2010-07-30
6 6
 Author: Rafael A Irizarry, Benilton S Carvalho <carvalho@bclab.org>, Robert Scharpf <rscharpf@jhsph.edu>, Matt Ritchie <mritchie@wehi.edu.au>
7 7
 Maintainer: Benilton S Carvalho <carvalho@bclab.org>, Robert Scharpf <rscharpf@jhsph.edu>, Matt Ritchie <mritchie@wehi.EDU.AU>
... ...
@@ -3111,7 +3111,51 @@ computeCopynumber.CNSet <- function(object, cnOptions){
3111 3111
 
3112 3112
 
3113 3113
 
3114
-
3115
-
3116
-
3117
-
3114
+## constructors for Illumina platform
3115
+constructIlluminaFeatureData <- function(gns, cdfName){
3116
+	pkgname <- paste(cdfName, "Crlmm", sep="")	
3117
+	path <- system.file("extdata", package=pkgname)
3118
+	load(file.path(path, "cnProbes.rda"))
3119
+	load(file.path(path, "snpProbes.rda"))
3120
+	cnProbes$chr <- chromosome2integer(cnProbes$chr)
3121
+	cnProbes <- as.matrix(cnProbes)
3122
+	snpProbes$chr <- chromosome2integer(snpProbes$chr)
3123
+	snpProbes <- as.matrix(snpProbes)
3124
+	mapping <- rbind(snpProbes, cnProbes, deparse.level=0)
3125
+	mapping <- mapping[match(gns, rownames(mapping)), ]
3126
+	isSnp <- 1L-as.integer(gns %in% rownames(cnProbes))
3127
+	mapping <- cbind(mapping, isSnp, deparse.level=0)
3128
+	stopifnot(identical(rownames(mapping), gns))
3129
+	colnames(mapping) <- c("chromosome", "position", "isSnp")
3130
+	new("AnnotatedDataFrame",
3131
+	    data=data.frame(mapping),
3132
+	    varMetadata=data.frame(labelDescription=colnames(mapping)))	
3133
+}
3134
+constructIlluminaAssayData <- function(np, snp, object, storage.mode="environment", order.index){
3135
+	stopifnot(identical(snp$gns, featureNames(object)))
3136
+	gns <- c(featureNames(object), np$gns)
3137
+	sns <- np$sns
3138
+	np <- np[1:2]
3139
+	snp <- snp[1:2]
3140
+	stripnames <- function(x) {
3141
+		dimnames(x) <- NULL
3142
+		x
3143
+	}
3144
+	np <- lapply(np, stripnames)
3145
+	snp <- lapply(snp, stripnames)
3146
+	A <- rbind(snp[[1]], np[[1]], deparse.level=0)[order.index, ]
3147
+	B <- rbind(snp[[2]], np[[2]], deparse.level=0)[order.index, ]
3148
+	gt <- stripnames(calls(object))
3149
+	emptyMatrix <- matrix(integer(), nrow(np[[1]]), ncol(A))
3150
+	gt <- rbind(gt, emptyMatrix, deparse.level=0)[order.index,]
3151
+	pr <- stripnames(snpCallProbability(object))
3152
+	pr <- rbind(pr, emptyMatrix, deparse.level=0)[order.index, ]
3153
+	emptyMatrix <- matrix(integer(), nrow(A), ncol(A))	
3154
+	aD <- assayDataNew(storage.mode,
3155
+			   alleleA=A,
3156
+			   alleleB=B,
3157
+			   call=gt,
3158
+			   callProbability=pr,
3159
+			   CA=emptyMatrix,
3160
+			   CB=emptyMatrix)
3161
+}
... ...
@@ -1,88 +0,0 @@
1
-setClassUnion("integerOrMissing", c("integer", "missing", "numeric"))
2
-setGeneric("totalCopyNumber", function(object, i, j, ...) standardGeneric("totalCopyNumber"))
3
-setMethod("totalCopyNumber",
4
-	  signature=signature(object="CNSet", i="integerOrMissing", j="integerOrMissing"),
5
-	  function(object, i, j, ...){
6
-	if(missing(i) & missing(j)){
7
-		if(inherits(CA(object), "ff") | inherits(CA(object), "ffdf")) stop("Must specify i and/or j for ff objects")
8
-	}
9
-	if(missing(i) & !missing(j)){
10
-		snp.index <- which(isSnp(object))	
11
-		cn.total <- as.matrix(CA(object)[, j])
12
-		if(length(snp.index) > 0){
13
-			cb <- as.matrix(CB(object)[snp.index, j])
14
-			snps <- (1:nrow(cn.total))[i %in% snp.index]
15
-			cn.total[snps, ] <- cn.total[snps, j] + cb				
16
-		}
17
-	}
18
-	if(!missing(i) & missing(j)){
19
-		snp.index <- intersect(which(isSnp(object)), i)
20
-		cn.total <- as.matrix(CA(object)[i, ])
21
-		if(length(snp.index) > 0){
22
-			cb <- as.matrix(CB(object)[snp.index, ])
23
-			snps <- (1:nrow(cn.total))[i %in% snp.index]
24
-			cn.total[snps, ] <- cn.total[snps, ] + cb				
25
-		}
26
-	}
27
-	if(!missing(i) & !missing(j)){
28
-		snp.index <- intersect(which(isSnp(object)), i)		
29
-		cn.total <- as.matrix(CA(object)[i, j])
30
-		if(length(snp.index) > 0){
31
-			cb <- as.matrix(CB(object)[snp.index, j])
32
-			snps <- (1:nrow(cn.total))[i %in% snp.index]
33
-			cn.total[snps, ] <- cn.total[snps, ] + cb
34
-		}
35
-	}
36
-	cn.total <- cn.total/100
37
-	dimnames(cn.total) <- NULL
38
-	return(cn.total)
39
-})
40
-
41
-
42
-constructFeatureData <- function(gns, cdfName){
43
-	pkgname <- paste(cdfName, "Crlmm", sep="")	
44
-	path <- system.file("extdata", package=pkgname)
45
-	load(file.path(path, "cnProbes.rda"))
46
-	load(file.path(path, "snpProbes.rda"))
47
-	cnProbes$chr <- chromosome2integer(cnProbes$chr)
48
-	cnProbes <- as.matrix(cnProbes)
49
-	snpProbes$chr <- chromosome2integer(snpProbes$chr)
50
-	snpProbes <- as.matrix(snpProbes)
51
-	mapping <- rbind(snpProbes, cnProbes, deparse.level=0)
52
-	mapping <- mapping[match(gns, rownames(mapping)), ]
53
-	isSnp <- 1L-as.integer(gns %in% rownames(cnProbes))
54
-	mapping <- cbind(mapping, isSnp, deparse.level=0)
55
-	stopifnot(identical(rownames(mapping), gns))
56
-	colnames(mapping) <- c("chromosome", "position", "isSnp")
57
-	new("AnnotatedDataFrame",
58
-	    data=data.frame(mapping),
59
-	    varMetadata=data.frame(labelDescription=colnames(mapping)))	
60
-}
61
-constructAssayData <- function(np, snp, object, storage.mode="environment", order.index){
62
-	stopifnot(identical(snp$gns, featureNames(object)))
63
-	gns <- c(featureNames(object), np$gns)
64
-	sns <- np$sns
65
-	np <- np[1:2]
66
-	snp <- snp[1:2]
67
-	stripnames <- function(x) {
68
-		dimnames(x) <- NULL
69
-		x
70
-	}
71
-	np <- lapply(np, stripnames)
72
-	snp <- lapply(snp, stripnames)
73
-	A <- rbind(snp[[1]], np[[1]], deparse.level=0)[order.index, ]
74
-	B <- rbind(snp[[2]], np[[2]], deparse.level=0)[order.index, ]
75
-	gt <- stripnames(calls(object))
76
-	emptyMatrix <- matrix(integer(), nrow(np[[1]]), ncol(A))
77
-	gt <- rbind(gt, emptyMatrix, deparse.level=0)[order.index,]
78
-	pr <- stripnames(snpCallProbability(object))
79
-	pr <- rbind(pr, emptyMatrix, deparse.level=0)[order.index, ]
80
-	emptyMatrix <- matrix(integer(), nrow(A), ncol(A))	
81
-	aD <- assayDataNew(storage.mode,
82
-			   alleleA=A,
83
-			   alleleB=B,
84
-			   call=gt,
85
-			   callProbability=pr,
86
-			   CA=emptyMatrix,
87
-			   CB=emptyMatrix)
88
-}
... ...
@@ -48,7 +48,9 @@ options(continue=" ")
48 48
 <<libraries>>=
49 49
 library(crlmm)
50 50
 crlmm:::validCdfNames()
51
-outdir <- "/thumper/ctsa/snpmicroarray/rs/data/hapmap/crlmmVignette/release/illumina"
51
+if(length(grep("development", sessionInfo()[[1]]$status)) == 1){
52
+	outdir <- "/thumper/ctsa/snpmicroarray/rs/data/hapmap/crlmmVignette/devel/illumina"	
53
+} else outdir <- "/thumper/ctsa/snpmicroarray/rs/data/hapmap/crlmmVignette/release/illumina"	
52 54
 dir.create(outdir, showWarnings=FALSE, recursive=TRUE)
53 55
 datadir <- "/thumper/ctsa/snpmicroarray/illumina/IDATS/370k"
54 56
 @ 
... ...
@@ -88,25 +90,24 @@ redfiles = all(file.exists(paste(arrayNames, "_Red.idat", sep="")))
88 90
 @ 
89 91
 
90 92
 <<samplesToProcess2, echo=FALSE>>=
91
-if(!exists("crlmmResult")){
92
-	if(!file.exists(file.path(outdir, "crlmmResult.rda"))){
93
-		RG <- readIdatFiles(samplesheet, 
94
-				    path=dirname(arrayNames[1]), 
95
-				    arrayInfoColNames=list(barcode=NULL, position="SentrixPosition"), 
96
-				    saveDate=TRUE)
97
-		crlmmResult <- crlmmIllumina(RG=RG, 
98
-					     cdfName="human370v1c", 
99
-					     sns=pData(RG)$ID, 
100
-					     returnParams=TRUE,
101
-					     cnFile=file.path(outdir, "cnFile.rda"),
102
-					     snpFile=file.path(outdir, "snpFile.rda"),
103
-					     save.it=TRUE)
104
-		protocolData(crlmmResult)$ScanDate <- protocolData(RG)$ScanDate
105
-		range(protocolData(crlmmResult)$ScanDate)
106
-		save(crlmmResult, file=file.path(outdir, "crlmmResult.rda"))
107
-		rm(RG); gc()
108
-	} 
109
-}
93
+## To speed up repeated calls to Sweave 
94
+RG <- checkExists("RG", .path=outdir,
95
+		  .FUN=readIdatFiles,
96
+		  sampleSheet=samplesheet,
97
+		  path=dirname(arrayNames[1]),
98
+		  arrayInfoColNames=list(barcode=NULL, position="SentrixPosition"), 
99
+		  saveDate=TRUE)
100
+annotation(RG) <- "human370v1c"
101
+crlmmResult <- checkExists("crlmmResult",
102
+			   .path=outdir,
103
+			   .FUN=crlmmIllumina,
104
+			   RG=RG,
105
+			   sns=pData(RG)$ID, 
106
+			   returnParams=TRUE,
107
+			   cnFile=file.path(outdir, "cnFile.rda"),
108
+			   snpFile=file.path(outdir, "snpFile.rda"),
109
+			   save.it=TRUE)
110
+protocolData(crlmmResult)$ScanDate <- protocolData(RG)$ScanDate
110 111
 @ 
111 112
 
112 113
 <<samplesToProcess3, eval=FALSE>>=
... ...
@@ -114,6 +115,7 @@ RG <- readIdatFiles(samplesheet,
114 115
 		    path=dirname(arrayNames[1]), 
115 116
 		    arrayInfoColNames=list(barcode=NULL, position="SentrixPosition"), 
116 117
 		    saveDate=TRUE)
118
+annotation(RG) <- "human370v1c"
117 119
 crlmmResult <- crlmmIllumina(RG=RG, 
118 120
 			     cdfName="human370v1c", 
119 121
 			     sns=pData(RG)$ID, 
... ...
@@ -122,8 +124,6 @@ crlmmResult <- crlmmIllumina(RG=RG,
122 124
 			     snpFile=file.path(outdir, "snpFile.rda"),
123 125
 			     save.it=TRUE)
124 126
 protocolData(crlmmResult)$ScanDate <- protocolData(RG)$ScanDate
125
-range(protocolData(crlmmResult)$ScanDate)
126
-rm(RG); gc()
127 127
 @ 
128 128
 
129 129
 \noindent Finally, we load a few of the intermediate files that were
... ...
@@ -133,17 +133,11 @@ load(file.path(outdir, "snpFile.rda"))
133 133
 res <- get("res")
134 134
 load(file.path(outdir, "cnFile.rda"))
135 135
 cnAB <- get("cnAB")
136
-load(file.path(outdir, "crlmmResult.rda"))
137 136
 @ 
138 137
 
139
-<<loadIntermediate, eval=TRUE, echo=FALSE>>=
140
-if(!exists("res")){
141
-	load(file.path(outdir, "snpFile.rda"))
142
-	res <- get("res")
143
-	load(file.path(outdir, "cnFile.rda"))
144
-	cnAB <- get("cnAB")
145
-	load(file.path(outdir, "crlmmResult.rda"))
146
-}
138
+<<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"))
147 141
 @ 
148 142
 
149 143
 After running the crlmm algorithm, we construct a container for