Browse code

Remove help file for ffdf-class. Adding documentation for CNSet accessors.

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

Rob Scharp authored on 04/10/2010 02:21:59
Showing7 changed files

... ...
@@ -1,13 +1,13 @@
1 1
 useDynLib("crlmm", .registration=TRUE)
2
-
3 2
 ## this is temporary
4 3
 ## exportPattern("^[^\\.]")
5 4
 
5
+##---------------------------------------------------------------------------
6 6
 ## Biobase
7
+##---------------------------------------------------------------------------
7 8
 importClassesFrom(Biobase, AnnotatedDataFrame, AssayData, eSet, SnpSet,
8 9
 		  NChannelSet, MIAME, Versioned, VersionedBiobase,
9 10
 		  Versions)
10
-
11 11
 importMethodsFrom(Biobase, annotation, "annotation<-",
12 12
                   annotatedDataFrameFrom, assayData, "assayData<-",
13 13
                   combine, dims, experimentData, "experimentData<-",
... ...
@@ -18,21 +18,19 @@ importMethodsFrom(Biobase, annotation, "annotation<-",
18 18
                   snpCallProbability,
19 19
 		  "snpCall<-", "snpCallProbability<-", storageMode,
20 20
                   "storageMode<-", updateObject, varLabels)
21
-
22 21
 importFrom(Biobase, assayDataElement, assayDataElementNames,
23 22
            assayDataElementReplace, assayDataNew, classVersion,
24 23
            validMsg)
25 24
 
25
+##---------------------------------------------------------------------------
26 26
 ## oligoClasses
27
+##---------------------------------------------------------------------------
27 28
 importClassesFrom(oligoClasses, SnpSuperSet, AlleleSet, CNSet)
28
-
29 29
 importMethodsFrom(oligoClasses, allele, calls, "calls<-", confs,
30 30
 		  "confs<-", cnConfidence, "cnConfidence<-", isSnp,
31 31
 		  chromosome, position, A, B,
32 32
 		  "A<-", "B<-", open, close, flags,
33 33
 		  batchStatistics, "batchStatistics<-", updateObject)
34
-
35
-
36 34
 importFrom(oligoClasses, chromosome2integer, celfileDate, list.celfiles,
37 35
            copyNumber, initializeBigMatrix, initializeBigVector)
38 36
 
... ...
@@ -41,21 +39,13 @@ importFrom(graphics, abline, axis, layout, legend, mtext, par, plot,
41 39
            polygon, rect, segments, text, points, boxplot, lines)
42 40
 
43 41
 importFrom(grDevices, grey)
44
-
45 42
 importFrom(affyio, read.celfile.header, read.celfile)
46
-
47 43
 importFrom(preprocessCore, normalize.quantiles.use.target, normalize.quantiles)
48
-
49 44
 importFrom(utils, data, packageDescription, setTxtProgressBar, txtProgressBar)
50
-
51 45
 importFrom(stats, coef, cov, dnorm, kmeans, lm, mad, median, quantile, sd, update)
52
-
53 46
 importFrom(genefilter, rowSds)
54
-
55 47
 importFrom(mvtnorm, dmvnorm)
56
-
57 48
 importFrom(ellipse, ellipse)
58
-
59 49
 importFrom(ff, ffdf, physical.ff, physical.ffdf, ffrowapply)
60 50
 ## It is important not to import these classes from oligoClasses
61 51
 ## Doing so causes the following errors:
... ...
@@ -64,7 +54,6 @@ importFrom(ff, ffdf, physical.ff, physical.ffdf, ffrowapply)
64 54
 ##  unable to find an inherited method for function "medianA.AA<-", for signature "CNSet", "ff_matrix"
65 55
 ##importClassesFrom(oligoClasses, ffdf, ff_matrix)
66 56
 importClassesFrom(oligoClasses, ff_matrix)
67
-
68 57
 ## Important to export these classes
69 58
 ##exportClasses(ff_or_matrix, ff_matrix, ffdf)
70 59
 ##exportClasses(ff_or_matrix)
... ...
@@ -73,7 +62,6 @@ exportMethods(CA, CB)
73 62
 export(crlmm,
74 63
        crlmmIllumina,
75 64
        crlmmIlluminaV2,
76
-##       ellipseCenters,
77 65
        genotype,
78 66
        readIdatFiles,
79 67
        snprma,
... ...
@@ -83,17 +71,7 @@ export(crlmm,
83 71
        crlmmCopynumber2, crlmmCopynumberLD, crlmmCopynumber)
84 72
 export(constructIlluminaCNSet)
85 73
 export(totalCopynumber)
86
-export(cnrma, cnrma2)
87 74
 exportMethods(A, B, nuA, nuB, phiA, phiB, corr, tau2, Ns, medians, mads)
88
-##export(genotypeSummary,
89
-##       shrinkSummary,
90
-##       estimateCnParameters,
91
-##       shrinkGenotypeSummaries,
92
-##       summarizeSnps,
93
-##       constructIlluminaAssayData,
94
-##       ACN, C1, C2, C3)
95
-## For debugging
96
-## exportPattern("^[^\\.]")
97 75
 
98 76
 
99 77
 
... ...
@@ -9,7 +9,7 @@ getProtocolData.Affy <- function(filenames){
9 9
 			                           row.names=colnames(scanDates)))
10 10
 	return(protocoldata)
11 11
 }
12
-getFeatureData.Affy <- function(cdfName, copynumber=FALSE){
12
+getFeatureData <- function(cdfName, copynumber=FALSE){
13 13
 	pkgname <- getCrlmmAnnotationName(cdfName)
14 14
 	if(!require(pkgname, character.only=TRUE)){
15 15
 		suggCall <- paste("library(", pkgname, ", lib.loc='/Altern/Lib/Loc')", sep="")
... ...
@@ -36,20 +36,22 @@ getFeatureData.Affy <- function(cdfName, copynumber=FALSE){
36 36
 	M <- matrix(NA, length(featurenames), 3, dimnames=list(featurenames, fvarlabels))
37 37
 	index <- match(rownames(snpProbes), rownames(M)) #only snp probes in M get assigned position
38 38
 	M[index, "position"] <- snpProbes[, grep("pos", colnames(snpProbes))]
39
-	M[index, "chromosome"] <- snpProbes[, grep("chr", colnames(snpProbes))]
39
+	M[index, "chromosome"] <- chromosome2integer(snpProbes[, grep("chr", colnames(snpProbes))])
40 40
 	M[index, "isSnp"] <- 1L
41 41
 	index <- which(is.na(M[, "isSnp"]))
42 42
 	M[index, "isSnp"] <- 1L
43 43
 	if(copynumber){
44 44
 		index <- match(rownames(cnProbes), rownames(M)) #only snp probes in M get assigned position
45 45
 		M[index, "position"] <- cnProbes[, grep("pos", colnames(cnProbes))]
46
-		M[index, "chromosome"] <- cnProbes[, grep("chr", colnames(cnProbes))]
46
+		M[index, "chromosome"] <- chromosome2integer(cnProbes[, grep("chr", colnames(cnProbes))])
47 47
 		M[index, "isSnp"] <- 0L
48 48
 	}
49 49
 	##A few of the snpProbes do not match -- I think it is chromosome Y.
50 50
 	M[is.na(M[, "isSnp"]), "isSnp"] <- 1L
51
-	return(new("AnnotatedDataFrame", data=data.frame(M)))
51
+	M <- data.frame(M)
52
+	return(new("AnnotatedDataFrame", data=M))
52 53
 }
54
+getFeatureData.Affy <- getFeatureData
53 55
 
54 56
 construct <- function(filenames,
55 57
 		      cdfName,
... ...
@@ -60,7 +62,7 @@ construct <- function(filenames,
60 62
 	}
61 63
 	if(missing(sns) & missing(filenames)) stop("one of filenames or samplenames (sns) must be provided")
62 64
 	if(verbose) message("Initializing container for copy number estimation")
63
-	featureData <- getFeatureData.Affy(cdfName, copynumber=copynumber)
65
+	featureData <- getFeatureData(cdfName, copynumber=copynumber)
64 66
 	if(!missing(fns)){
65 67
 		index <- match(fns, featureNames(featureData))
66 68
 		if(all(is.na(index))) stop("fns not in featureNames")
... ...
@@ -118,7 +120,7 @@ genotype <- function(filenames,
118 120
 			     verbose=verbose,
119 121
 			     batch=batch)
120 122
 	##save(callSet, file=file.path(outdir, "callSet.rda"))
121
-	open(callSet)
123
+	if(is.lds) open(callSet)
122 124
 	mixtureParams <- matrix(NA, 4, length(filenames))
123 125
 	is.snp <- isSnp(callSet)
124 126
 	snp.index <- which(is.snp)
... ...
@@ -160,6 +162,7 @@ genotype <- function(filenames,
160 162
 		       cnrma2=cnrma2(...))
161 163
 	}
162 164
 	## consider passing only A for NPs.
165
+	if(verbose) message("Quantile normalizing nonpolymorphic markers")
163 166
 	AA <- cnrmaFxn(FUN, A=A(callSet),
164 167
 		       filenames=filenames,
165 168
 		       row.names=featureNames(callSet)[np.index],
... ...
@@ -167,7 +170,10 @@ genotype <- function(filenames,
167 170
 		       sns=sns,
168 171
 		       seed=seed,
169 172
 		       verbose=verbose)
173
+	if(verbose) message("Saving callSet.rda")
174
+	save(callSet, file=file.path(outdir, "callSet.rda"))
170 175
 	if(!is.lds) A(callSet) <- AA
176
+	## otherwise, the normalized values were written to file... no need to do anything
171 177
 	rm(AA)
172 178
 	FUN <- ifelse(is.lds, "crlmmGT2", "crlmmGT")
173 179
 	## genotyping
... ...
@@ -176,6 +182,7 @@ genotype <- function(filenames,
176 182
 		       crlmmGT2=crlmmGT2(...),
177 183
 		       crlmmGT=crlmmGT(...))
178 184
 	}
185
+	if(verbose) message("Call genotypes at polymorphic markers")
179 186
 	tmp <- crlmmGTfxn(FUN,
180 187
 			  A=snprmaRes[["A"]],
181 188
 			  B=snprmaRes[["B"]],
... ...
@@ -1732,8 +1739,6 @@ summarizeSnps <- function(strata,
1732 1739
 	if(is.lds) return(TRUE) else return(object)
1733 1740
 }
1734 1741
 
1735
-
1736
-
1737 1742
 crlmmCopynumber <- function(object,
1738 1743
 			    MIN.SAMPLES=10,
1739 1744
 			    SNRMin=5,
... ...
@@ -1,5 +1,5 @@
1 1
 ## Need to submit from a fat node!
2
-all:	illumina_copynumber copynumber 
2
+all:	illumina affy 
3 3
 
4 4
 cleanScripts:
5 5
 	-rm *.pdf *.tex *.log *.aux 
... ...
@@ -24,19 +24,16 @@ testingFF: testingFF.Rnw
24 24
 	cat ~/bin/cluster.template | perl -pe "s/Rprog/testingFF.R/" > testingFF.R.sh
25 25
 	qsub -m e -r y -cwd -l mem_free=12G,h_vmem=16G testingFF.R.sh
26 26
 
27
-copynumber:	copynumber.Rnw
27
+affy:	copynumber.Rnw
28 28
 	echo "Stangle(\"copynumber.Rnw\")" | R --no-save --no-restore;	
29 29
 	cat ~/bin/cluster.template | perl -pe "s/Rprog/copynumber.R/" > copynumber.R.sh
30 30
 	qsub -m e -r y -cwd -l mem_free=12G,h_vmem=16G copynumber.R.sh
31 31
 
32
-illumina_copynumber:	illumina_copynumber.Rnw
32
+illumina:	illumina_copynumber.Rnw
33 33
 	echo "Stangle(\"illumina_copynumber.Rnw\")" | R --no-save --no-restore;	
34 34
 	cat ~/bin/cluster.template | perl -pe "s/Rprog/illumina_copynumber.R/" > illumina_copynumber.R.sh
35 35
 	qsub -m e -r y -cwd -l mem_free=12G,h_vmem=16G illumina_copynumber.R.sh
36 36
 
37
-sweavecopynumber: copynumber.Rnw
38
-	echo "Sweave(\"copynumber.Rnw\")" | R --no-save --no-restore;	
39
-
40 37
 texcopynumber:	copynumber.tex
41 38
 	texi2dvi --pdf --clean --quiet copynumber.tex
42 39
 
... ...
@@ -56,9 +56,9 @@ options(continue=" ")
56 56
 @
57 57
 
58 58
 <<libraries>>=
59
-##library(ff)
60
-library(crlmm)
61
-crlmm:::validCdfNames()
59
+library(ff)
60
+pkgs <- annotationPackages()
61
+pkgs[grep("Crlmm", pkgs)]
62 62
 if(getRversion() < "2.12.0"){
63 63
 	rpath <- getRversion()
64 64
 } else rpath <- "trunk"
... ...
@@ -84,7 +84,7 @@ these directories exist.
84 84
 
85 85
 <<checkSetup>>=
86 86
 if(!file.exists(outdir)) stop("Please specify valid directory for storing output")
87
-if(!file.exists(pathToCels)) stop("Please specify the correct path to the CEL files")
87
+if(!file.exists(datadir)) stop("Please specify the correct path to the CEL files")
88 88
 @
89 89
 
90 90
 \paragraph{Preprocessing Illumina IDAT files.}
... ...
@@ -123,110 +123,158 @@ grnfiles = all(file.exists(paste(arrayNames, "_Grn.idat", sep="")))
123 123
 redfiles = all(file.exists(paste(arrayNames, "_Red.idat", sep="")))
124 124
 @
125 125
 
126
-<<eval=FALSE>>=
127
-RG <- checkExists("RG", .path=outdir,
128
-		  .FUN=readIdatFiles,
129
-		  sampleSheet=samplesheet,
130
-		  path=dirname(arrayNames[1]),
131
-		  arrayInfoColNames=list(barcode=NULL, position="SentrixPosition"),
132
-		  saveDate=TRUE)
133
-annotation(RG) <- "human370v1c"
134
-crlmmResult <- checkExists("crlmmResult",
135
-			   .path=outdir,
136
-			   .FUN=crlmmIllumina,
137
-			   RG=RG,
138
-			   sns=pData(RG)$ID,
139
-			   returnParams=TRUE,
140
-			   cnFile=file.path(outdir, "cnFile.rda"),
141
-			   snpFile=file.path(outdir, "snpFile.rda"),
142
-			   save.it=TRUE)
143
-protocolData(crlmmResult)$ScanDate <- protocolData(RG)$ScanDate
126
+<<genotype>>=
127
+load.it <- TRUE
128
+if(!load.it){
129
+	rmFiles <- list.files(outdir, pattern=".ff", full.names=TRUE)
130
+	unlink(rmFiles)
131
+}
132
+container <- checkExists("container", .path=outdir,
133
+			 .FUN=crlmm:::genotype.Illumina,
134
+			 sampleSheet=samplesheet,
135
+			 path=dirname(arrayNames[1]),
136
+			 arrayInfoColNames=list(barcode=NULL, position="SentrixPosition"),
137
+			 cdfName="human370v1c",
138
+			 .load.it=load.it)
139
+cnSet <- checkExists("cnSet", .path=outdir, .FUN=crlmmCopynumber, object=container, .load.it=load.it)
144 140
 @
145 141
 
146
-<<samplesToProcess2, echo=FALSE>>=
147
-crlmmResult <- checkExists("crlmmResult",
148
-			   .path=outdir,
149
-			   .FUN=crlmmIlluminaV2,
150
-			   sampleSheet=samplesheet,
151
-			   path=dirname(arrayNames[1]),
152
-			   arrayInfoColNames=list(barcode=NULL, position="SentrixPosition"),
153
-			   returnParams=TRUE,
154
-			   cnFile=file.path(outdir, "cnFile.rda"),
155
-			   snpFile=file.path(outdir, "snpFile.rda"),
156
-			   save.it=TRUE,
157
-			   cdfName="human370v1c", .load.it=FALSE)
158
-## MR: shouldn't this be intensities in the A element of res?
159
-load(file.path(outdir, "snpFile.rda"))
160
-res[["A"]][1:5, 1:5]
161
-stop()
142
+<<copynumberObject>>=
143
+marker.index <- which(chromosome(cnSet) <= 22 & isSnp(cnSet))
144
+sample.index <- 1:5 ## first five samples
145
+invisible(open(cnSet))
146
+copynumberSet <- as(cnSet[marker.index, 1:5], "CopyNumberSet")
147
+invisible(close(cnSet))
148
+copynumberSet <- copynumberSet[order(chromosome(copynumberSet), position(copynumberSet)), ]
149
+indices <- split(1:nrow(copynumberSet), chromosome(copynumberSet))
150
+dup.index <- unlist(sapply(indices, function(i, position)  i[duplicated(position[i])], position=position(copynumberSet)))
151
+if(length(dup.index) > 0) copynumberSet <- copynumberSet[-dup.index, ]
152
+## exclude any rows that are all missing
153
+missing.index <- which(rowSums(is.na(copyNumber(copynumberSet))) == ncol(copynumberSet))
154
+copynumberSet <- copynumberSet[-missing.index, ]
162 155
 @
163 156
 
164
-%\noindent Finally, we load a few of the intermediate files that were
165
-%created during the preprocessing and genotyping.
166
-%<<loadIntermediate, eval=FALSE>>=
167
-%load(file.path(outdir, "snpFile.rda"))
168
-%res <- get("res")
169
-%load(file.path(outdir, "cnFile.rda"))
170
-%cnAB <- get("cnAB")
171
-%@
172
-%
173
-%<<loadIntermediate, echo=FALSE>>=
174
-%trace(checkExists, browser)
175
-%res <- checkExists("snpFile", .path=outdir, .FUN=load, file=file.path(outdir, "snpFile.rda"))
176
-%cnAB <- checkExists("cnFile", .path=outdir, .FUN=load, file=file.path(outdir, "cnFile.rda"))
177
-%@
178
-
179
-After running the crlmm algorithm, we construct a container for
180
-storing the quantile normalized intensities, genotype calls, and
181
-allele-specific copy number estimates.  A few helper functions for
182
-facilitating the construction of this container have been added to the
183
-inst/scripts directory of this package and can be sourced as follows.
184
-
185
-<<constructContainer2>>=
186
-load(file.path(outdir, "snpFile.rda"))
187
-batch <- as.factor(rep("1", ncol(crlmmResult)))
188
-container <- constructIlluminaCNSet(crlmmResult=crlmmResult, snpFile=file.path(outdir, "snpFile.rda"), cnFile=file.path(outdir, "cnFile.rda"), batch=batch)
189
-A(container)[1:5, 1:5]
157
+<<hmm>>=
158
+if(require(VanillaICE)){
159
+	hmmOpts <- hmm.setup(copynumberSet, c("hom-del", "hem-del", "normal", "amp"),
160
+			     copynumberStates=0:3, normalIndex=3,
161
+			     log.initialP=rep(log(1/4), 4))
162
+	timing <- system.time(fit.cn <- hmm(copynumberSet, hmmOpts, verbose=FALSE))
163
+	hmm.df <- as.data.frame(fit.cn)
164
+	print(hmm.df[1:5, c(2:4,7:9)])
165
+}
190 166
 @
191 167
 
192
-As described in the \texttt{copynumber} vignette, two \R{} functions
193
-for copy number estimation are available: \Rfunction{crlmmCopynumber}
194
-and \Rfunction{crlmmCopynumber2}.  The latter requires that the assay
195
-data elements are represented using \Robject{ff} objects.  As the
196
-dataset for this vignette is small (43 arrays) and the above steps did
197
-not make use of the \Robject{ff} features, constructing \Robject{ff}
198
-objects at this point in the analysis would serve little purpose.  The
199
-decision to use ordinary matrices or \Robject{ff} objects should be
200
-decided at the beginning of the analysis and then propogated to both
201
-the genotyping and copy number estimation steps.  Here, we use the
202
-\Rfunction{crlmmCopynumber} to estimate copy number.
203
-
204
-<<crlmmCopynumber>>=
205
-##trace(indexComplete, browser)
206
-##marker.index <- which(isSnp(container) & chromosome(container) < 23)
207
-##shrinkGenotypeSummaries(container, index.list=list(marker.index), strata=1, is.lds=T, verbose=T, MIN.OBS=1, DF.PRIOR=50, MIN.SAMPLES=10)
208
-##summarizeSnps(strata=1, index.list=list(marker.index), object=container, GT.CONF.THR=0.95, verbose=T, is.lds=T, CHR.X=F)
209
-cnSet <- checkExists("cnSet", .path=outdir, .FUN=crlmmCopynumber, object=container)
168
+<<cbs>>=
169
+library("DNAcopy")
170
+CNA.object <- CNA(genomdat=copyNumber(copynumberSet),
171
+		  chrom=chromosome(copynumberSet),
172
+		  maploc=position(copynumberSet),
173
+		  data.type="logratio",
174
+		  sampleid=sampleNames(copynumberSet))
175
+smu.object <- smooth.CNA(CNA.object)
176
+if(!load.it) unlink(file.path(outdir, "cbs.segments.rda"))
177
+cbs.segments <- checkExists("cbs.segments", .path=outdir, .FUN=segment, x=smu.object, .load.it=load.it)
178
+cbs.segments <- cbind(cbs.segments$output, cbs.segments$segRows)
179
+cbs.segments$call <- rep(3, nrow(cbs.segments))
180
+cbs.segments$call[cbs.segments$seg.mean > 2.5] <- 4
181
+cbs.segments$call[cbs.segments$seg.mean < 1.25 & cbs.segments$seg.mean > 0.75] <- 2
182
+cbs.segments$call[cbs.segments$seg.mean < 0.75] <- 1
183
+cbs.ir <- RangedData(IRanges(cbs.segments$loc.start, cbs.segments$loc.end),
184
+		     chrom=cbs.segments$chrom,
185
+		     numMarkers=cbs.segments$num.mark,
186
+		     seg.mean=cbs.segments$seg.mean,
187
+		     cnCall=cbs.segments$call,
188
+		     sampleId=cbs.segments$ID)
210 189
 @
211 190
 
212
-<<estimateCopynumber, eval=FALSE>>=
213
-cnSet <- crlmmCopynumberLD(container)
191
+Plot the total copy number versus physical position for chromosome 1.
192
+Overlay the copy number predictions from the HMM and a CNV call from
193
+the segmentation obtained by thresholding the segment means.
194
+
195
+<<smoothFig, fig=TRUE, width=8, height=8>>=
196
+require(SNPchip)
197
+CHR <- 1
198
+cnSet2 <- copynumberSet
199
+fit.cn2 <- fit.cn
200
+cbs.ir2 <- cbs.ir
201
+copynumberSet <- copynumberSet[which(chromosome(copynumberSet) == CHR), ]
202
+layout(matrix(1:5, 5,1), heights=c(1, 1, 1, 1, 1))
203
+par(mar=c(0.5, 0.5, 0.5, 0.5), las=1, oma=c(4, 4, 4, 4))
204
+for(j in 1:5){
205
+	fit.cn <- fit.cn2[fit.cn2$chrom==CHR & fit.cn2$sampleId == sampleNames(copynumberSet)[j], ]
206
+	cbs.ir <- cbs.ir2[cbs.ir2$chrom==CHR & cbs.ir2$sampleId == sampleNames(copynumberSet)[j], ]
207
+	plot(position(copynumberSet), copyNumber(copynumberSet)[, j], pch=".", cex=0.6,col="grey60", xaxt="n",
208
+	     ylim=c(-1,6), ylab="total copy number",
209
+	     yaxt="n",
210
+	     xlim=c(0, max(position(copynumberSet))))
211
+	axis(2, at=0:6, labels=0:6)
212
+	if(j == 5){
213
+		mtext("physical position (Mb)", 1, outer=TRUE, line=2)
214
+		mtext("Chromosome 8", 3, line=0, outer=TRUE)
215
+		at <- pretty(position(copynumberSet))
216
+	}
217
+	##Use polygon at bottom to show predictions
218
+	w <- width(fit.cn)
219
+	fit.cn <- fit.cn[order(w, decreasing=TRUE), ]
220
+	##fit.cn <- fit.cn[fit.cn$numMarkers > 5, ]
221
+	col <- c("blue", "black", "grey80", "red")
222
+	x <- c(min(start(fit.cn)), max(end(fit.cn)))
223
+	xx <- c(x, rev(x))
224
+	y <- rep(c(-0.3, 0), each=2)
225
+	polygon(xx, y, col="white")
226
+	for(i in seq_along(w)){
227
+		x <- c(start(fit.cn)[i], end(fit.cn)[i])
228
+		xx <- c(x, rev(x))
229
+		statecolor <- col[fit.cn$state[i]]
230
+		polygon(xx, y, col=statecolor, border=statecolor)
231
+	}
232
+	w <- width(cbs.ir)
233
+	cbs.ir <- cbs.ir[order(w, decreasing=TRUE), ]
234
+	##cbs.ir <- cbs.ir[cbs.ir$numMarkers > 5, ]
235
+	y <- rep(c(-1, -0.7), each=2)
236
+	for(i in seq_along(w)){
237
+		x <- c(start(cbs.ir)[i], end(cbs.ir)[i])
238
+		xx <- c(x, rev(x))
239
+		statecolor <- col[cbs.ir$cnCall[i]]
240
+		polygon(xx, y, col=statecolor, border=statecolor)
241
+	}
242
+	axis(4, at=c(-0.9, -0.1), c("CBS", "HMM"), cex=0.8)
243
+	data("chromosomeAnnotation")
244
+	x <- as.integer(chromosomeAnnotation[8, 1:2])
245
+	x <- c(x, rev(x))
246
+	y <- c(-1, -1, 6, 6)
247
+	polygon(x, y, border="white", col="white")
248
+	## draw centromere
249
+	if(j==5){
250
+		legend("topright", fill=col, legend=c(0, 1, 2, expression(3+.)))
251
+##		invisible(plotCytoband(8, label.cytoband=FALSE, xlim=c(0,max(position(copynumberSet)))))
252
+		labels <- at/1e6
253
+		axis(1, at=at, labels=labels, outer=T, line=0)
254
+	}
255
+}
214 256
 @
215 257
 
216
-%In the following code, we define helper functions to construct a
217
-%\Robject{featureData} object that chromosome and physical position
218
-%annotation for each marker. We then define a constructor for
219
-%initializing the assay data -- the copy number estimation algorithm
220
-%requires normalized intensities, genotype calls, and genotype
221
-%confidence scores.  Note that the order of the construction is
222
-%important.  In particular, the validity method for the
223
-%\Robject{CNSetLM} object requires a 'batch' label in the
224
-%\Robject{protocolData} slot and that 'chromosome', 'position', and
225
-%'isSnp' labels are present in the \Robject{featureData} slot.  With
226
-%the object \Robject{cnSet} in hand, one can proceed with the copy
227
-%number analysis as outlined in the \texttt{copynumber.pdf} vignette.
228
-%Useful accessors and visualizations of the locus-level estimates are
229
-%also discussed in the \texttt{copynumber.pdf} vignette.
258
+
259
+Make a 4 x 4 table of the called states for each sample.
260
+
261
+<<smooth>>=
262
+fit.cn <- fit.cn2; cbs.ir <- cbs.ir2
263
+nm.hmm <- sapply(split(fit.cn$numMarkers, fit.cn$sampleId), sum)
264
+stopifnot(length(unique(nm.hmm)) == 1)
265
+nm.cbs <- sapply(split(cbs.ir$numMarkers, cbs.ir$sampleId), sum)
266
+stopifnot(length(unique(nm.cbs)) == 1)
267
+stopifnot(all.equal(nm.cbs, nm.hmm))
268
+hmm.states <- rep(fit.cn$state, fit.cn$numMarkers)
269
+sample.id <- rep(fit.cn$sampleId, fit.cn$numMarkers)
270
+hmm.states <- split(hmm.states, sample.id)
271
+cbs.states <- rep(cbs.ir$cnCall, cbs.ir$numMarkers)
272
+sample.id <- rep(cbs.ir$sampleId, cbs.ir$numMarkers)
273
+cbs.states <- split(cbs.states, sample.id)
274
+tabs <- vector("list", length(cbs.states))
275
+for(i in seq_along(cbs.states)) tabs[[i]] <- table(cbs.states[[i]], hmm.states[[i]])
276
+names(tabs) <- names(cbs.states)
277
+@
230 278
 
231 279
 
232 280
 \paragraph{Accessors for extracting the locus-level copy number
... ...
@@ -238,8 +286,8 @@ number for polymorphic markers on chromosome 21.
238 286
 
239 287
 <<acn-accessor>>=
240 288
 marker.index <- which(chromosome(cnSet) == 21 & isSnp(cnSet))
241
-ca <- CA(cnSet)[marker.index, ]/100
242
-cb <- CB(cnSet)[marker.index, ]/100
289
+ca <- CA(cnSet, i=marker.index)
290
+cb <- CB(cnSet, i=marker.index)
243 291
 missing.index <- which(rowSums(is.na(ca))==ncol(cnSet))
244 292
 ca <- ca[-missing.index, ]
245 293
 cb <- cb[-missing.index, ]
... ...
@@ -247,7 +295,7 @@ cb <- cb[-missing.index, ]
247 295
 \noindent Negating the \Robject{isSnp} function could be used to
248 296
 extract the estimates at nonpolymorphic markers. For instance,
249 297
 <<monomorphic-accessor>>=
250
-cn.monomorphic <- CA(cnSet)[which(chromosome(cnSet) == 21 & !isSnp(cnSet)), ]/100
298
+cn.monomorphic <- CA(cnSet, i=which(chromosome(cnSet) == 21 & !isSnp(cnSet)))
251 299
 @
252 300
 
253 301
 At polymorphic loci, the total copy number is the sum of the number of
... ...
@@ -260,9 +308,8 @@ available in the devel version of the \Rpackage{oligoClasses}.
260 308
 
261 309
 <<copyNumberHelper>>=
262 310
 cn.total <- ca+cb
263
-cn.total2 <- totalCopyNumber(cnSet, i=marker.index)
311
+cn.total2 <- totalCopynumber(cnSet, i=marker.index)
264 312
 cn.total2 <- cn.total2[-missing.index, ]
265
-dimnames(cn.total) <- NULL
266 313
 all.equal(cn.total2, cn.total)
267 314
 @
268 315
 
... ...
@@ -276,12 +323,14 @@ samples with the lowest (top) and highest (bottom) signal to noise
276 323
 ratios.
277 324
 
278 325
 <<snr, fig=TRUE>>=
279
-hist(cnSet$SNR, breaks=15)
326
+open(cnSet$SNR)
327
+snr <- cnSet$SNR[,]
328
+hist(snr, breaks=15)
280 329
 @
281 330
 
282 331
 <<firstSample, fig=TRUE>>=
283
-low.snr <- which(cnSet$SNR == min(cnSet$SNR))
284
-high.snr <- which(cnSet$SNR == max(cnSet$SNR))
332
+low.snr <- which(snr == min(snr))
333
+high.snr <- which(snr == max(snr))
285 334
 x <- position(cnSet)[marker.index]
286 335
 x <- x[-missing.index]
287 336
 par(mfrow=c(2,1), las=1, mar=c(0.5, 4, 0.5, 0.5), oma=c(4, 1, 1,1))
... ...
@@ -5,10 +5,14 @@
5 5
 \alias{CB,CNSet-method}
6 6
 \alias{lines,CNSet-method}
7 7
 \alias{totalCopynumber,CNSet-method}
8
+\alias{nuA,CNSet-method}
9
+\alias{nuB,CNSet-method}
8 10
 \alias{Ns,CNSet-method}
9 11
 \alias{corr,CNSet-method}
10 12
 \alias{mads,CNSet-method}
11 13
 \alias{medians,CNSet-method}
14
+\alias{phiA,CNSet-method}
15
+\alias{phiB,CNSet-method}
12 16
 \alias{tau2,CNSet-method}
13 17
 
14 18
 \title{crlmm methods for class "CNSet"}
... ...
@@ -32,6 +36,10 @@
32 36
     \item{CB}{\code{signature(object="CNSet")}: ...}
33 37
     \item{lines}{\code{signature(x="CNSet")}: ...}
34 38
     \item{totalCopynumber}{\code{signature(object="CNSet")}: ...}
39
+    \item{nuA}{\code{signature(object="CNSet")}: ...}
40
+    \item{nuB}{\code{signature(object="CNSet")}: ...}
41
+    \item{phiA}{\code{signature(object="CNSet")}: ...}
42
+    \item{phiB}{\code{signature(object="CNSet")}: ...}
35 43
     \item{Ns}{\code{signature(object="CNSet")}: ...}
36 44
     \item{corr}{\code{signature(object="CNSet")}: ...}
37 45
     \item{mads}{\code{signature(x="CNSet")}: ...}
... ...
@@ -1,6 +1,10 @@
1 1
 \name{copynumberAccessors}
2 2
 \alias{CA}
3 3
 \alias{CB}
4
+\alias{nuA}
5
+\alias{nuB}
6
+\alias{phiA}
7
+\alias{phiB}
4 8
 \alias{totalCopynumber}
5 9
 
6 10
 \title{
... ...
@@ -51,8 +55,10 @@ totalCopynumber(object,...)
51 55
 }
52 56
 
53 57
 \examples{
58
+## Version 1.6* of crlmm used CNSetLM objects.
54 59
 data(sample.CNSetLM)
55
-## update to class CNSet
60
+
61
+## To update to class CNSet, use
56 62
 cnSet <- as(sample.CNSetLM, "CNSet")
57 63
 all(isCurrent(cnSet)) ## is the cnSet object current?
58 64
 
59 65
deleted file mode 100644
... ...
@@ -1,31 +0,0 @@
1
-\name{ffdf-class}
2
-\Rdversion{1.1}
3
-\docType{class}
4
-\alias{ffdf-class}
5
-\alias{annotatedDataFrameFrom,ff_matrix-method}
6
-\alias{list-class}
7
-\title{Class "ffdf"}
8
-\description{
9
-	object of class \code{ffdf}
10
-}
11
-\section{Objects from the Class}{A virtual Class: No objects may be created from it.}
12
-\section{Slots}{
13
-  \describe{
14
-    \item{\code{.S3Class}:}{Object of class \code{"character"} ~~ }
15
-  }
16
-}
17
-\section{Extends}{
18
-Class \code{"\linkS4class{oldClass}"}, directly.
19
-Class \code{"\linkS4class{list_or_ffdf}"}, directly.
20
-Class \code{"\linkS4class{ff_or_matrix}"}, directly.
21
-}
22
-\section{Methods}{
23
-No methods defined with class "ffdf" in the signature.
24
-}
25
-\seealso{
26
-	\code{link{ffdf}}
27
-}
28
-\examples{
29
-showClass("ffdf")
30
-}
31
-\keyword{classes}