Browse code

updated first half of inst/scripts/copynumber.Rnw, as well as a few of the help files

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

Rob Scharp authored on 23/12/2009 10:57:10
Showing16 changed files

... ...
@@ -350,3 +350,7 @@ is decoded and scanned
350 350
 2009-12-04 B. Carvalho - committed version 1.5.11
351 351
 
352 352
  * moved IRanges and oligoClasses to Depends (DESCRIPTION)
353
+
354
+2009-12-23 R. Scharpf - commited version 1.5.17
355
+
356
+** updated first half of inst/scripts/copynumber.Rnw, as well as a few of the help files
... ...
@@ -1,14 +1,14 @@
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.5.16
5
-Date: 2009-12-03
4
+Version: 1.5.17
5
+Date: 2009-12-23
6 6
 Author: Rafael A Irizarry, Benilton S Carvalho <bcarvalh@jhsph.edu>, Robert Scharpf <rscharpf@jhsph.edu>, Matt Ritchie <mritchie@wehi.edu.au>
7 7
 Maintainer: Benilton S Carvalho <bcarvalh@jhsph.edu>, Robert Scharpf <rscharpf@jhsph.edu>, Matt Ritchie <mritchie@wehi.EDU.AU>
8 8
 Description: Faster implementation of CRLMM specific to SNP 5.0 and 6.0 arrays, as well as a copy number tool specific to 5.0, 6.0, and Illumina platforms
9 9
 License: Artistic-2.0
10
-Depends: methods, Biobase (>= 2.7.2), R (>= 2.11.0)
11
-Imports: affyio, preprocessCore, utils, stats, genefilter, splines, mvtnorm, oligoClasses (>= 1.9.20), ellipse, methods, SNPchip, IRanges (>= 1.5.21)
10
+Depends: methods, Biobase (>= 2.7.2), R (>= 2.11.0), oligoClasses (>= 1.9.21)
11
+Imports: affyio, preprocessCore, utils, stats, genefilter, splines, mvtnorm, ellipse, SNPchip
12 12
 Suggests: hapmapsnp5, hapmapsnp6, genomewidesnp5Crlmm (>= 1.0.2),genomewidesnp6Crlmm (>= 1.0.2), snpMatrix
13 13
 Collate: AllGenerics.R
14 14
          methods-AlleleSet.R
... ...
@@ -22,13 +22,14 @@ importFrom(Biobase, assayDataElement, assayDataElementNames,
22 22
 importClassesFrom(oligoClasses, SnpSuperSet, AlleleSet)
23 23
 
24 24
 importMethodsFrom(oligoClasses, allele, calls, "calls<-", confs,
25
-		  "confs<-", cnConfidence, "cnConfidence<-", copyNumber,
26
-		  isSnp, chromosome, position, CA, "CA<-", CB, "CB<-",
27
-		  rangedData, segmentData, "rangedData<-",
28
-		  "segmentData<-")
25
+		  "confs<-", cnConfidence, "cnConfidence<-", 
26
+		  isSnp, chromosome, position, CA, "CA<-", CB, "CB<-")
27
+
29 28
 
30 29
 importFrom(oligoClasses, chromosome2integer, celfileDate, list.celfiles)
31 30
 
31
+importFrom(oligoClasses, copyNumber)
32
+
32 33
 importFrom(graphics, abline, axis, layout, legend, mtext, par, plot,
33 34
            polygon, rect, segments, text, points, boxplot)
34 35
 
... ...
@@ -48,18 +49,13 @@ importFrom(mvtnorm, dmvnorm)
48 49
 
49 50
 importFrom(ellipse, ellipse)
50 51
 
51
-## move to oligoClasses
52
-export(celDates, crlmm, cnOptions, calls, "calls<-", snprma)
53
-
54
-##exportMethods(computeHmm) ##move to VanillaICE
55 52
 exportMethods(A, B, copyNumber)
56
-export(crlmmCopynumber) 
57
-export(ellipse) 
58
-export(##viterbi.CNSet, ##move to VanillaICE
59
-	combineIntensities,
60
-	whichPlatform,
61
-	isValidCdfName,
62
-	crlmmWrapper,
63
-	readIdatFiles,
64
-        withinGenotypeMoments,  
65
-	locationAndScale, getParam.SnpSuperSet)
53
+export(cnOptions, crlmm, crlmmCopynumber, ellipse, readIdatFiles, snprma) 
54
+
55
+##export(##viterbi.CNSet, ##move to VanillaICE
56
+##	combineIntensities,
57
+##	whichPlatform,
58
+##	isValidCdfName,
59
+##	crlmmWrapper,
60
+##        withinGenotypeMoments,  
61
+##	locationAndScale, getParam.SnpSuperSet)
... ...
@@ -15,3 +15,4 @@ setGeneric("pr", function(object, name, batch, value) standardGeneric("pr"))
15 15
 setGeneric("snpIndex", function(object) standardGeneric("snpIndex"))
16 16
 setGeneric("snpNames", function(object) standardGeneric("snpNames"))
17 17
 ##setGeneric("splitByChromosome", function(object, ...) standardGeneric("splitByChromosome"))
18
+
... ...
@@ -225,7 +225,9 @@ crlmmCopynumber <- function(filenames, cnOptions, object, ...){
225 225
 		k <- grep("chr", colnames(snpProbes))
226 226
 		if(length(k) < 1) stop("chr or chromosome not in colnames(snpProbes)")
227 227
 	}
228
+	set.seed(cnOptions[["seed"]])  ##for reproducibility
228 229
 	chromosome <- cnOptions[["chromosome"]]
230
+	SNRmin <- cnOptions[["SNRmin"]]
229 231
 	for(CHR in chromosome){
230 232
 		cat("Chromosome ", CHR, "\n")
231 233
 		if(createIntermediateObjects){
... ...
@@ -247,6 +249,10 @@ crlmmCopynumber <- function(filenames, cnOptions, object, ...){
247 249
 			cnSet <- combineIntensities(res=snpI,
248 250
 						    NP=npI,
249 251
 						    callSet=callSet[index.snps, ])
252
+			if(any(cnSet$SNR > SNRmin)){
253
+				warning(paste("Excluding samples with SNR < ", SNRmin))
254
+				cnSet <- cnSet[, cnSet$SNR >= SNRmin]
255
+			}
250 256
 			rm(snpI, npI, snps, cnps, index.snps, index.nps); gc()
251 257
 			pData(cnSet)$batch <- cnOptions[["batch"]]
252 258
 			featureData(cnSet) <- lm.parameters(cnSet, cnOptions)
... ...
@@ -262,10 +268,13 @@ crlmmCopynumber <- function(filenames, cnOptions, object, ...){
262 268
 			rm(cnSet, alleleSet); gc()
263 269
 			next()
264 270
 		}
265
-		if(cnOptions[["save.cnset"]]){
266
-			save(cnSet, file=file.path(cnOptions[["outdir"]], paste("cnSet_", CHR, ".rda", sep="")))
271
+		if(length(chromosome) == 1){
272
+			if(cnOptions[["save.cnset"]]){
273
+				save(cnSet, file=file.path(cnOptions[["outdir"]], paste("cnSet_", CHR, ".rda", sep="")))
274
+			}
267 275
 		}
268 276
 		if(length(chromosome) > 1){
277
+			save(cnSet, file=file.path(cnOptions[["outdir"]], paste("cnSet_", CHR, ".rda", sep="")))
269 278
 			rm(cnSet); gc()
270 279
 		} else {
271 280
 			return(cnSet)
... ...
@@ -287,7 +296,6 @@ crlmmWrapper <- function(filenames, cnOptions, ...){
287 296
 	rgFile=cnOptions[["rgFile"]]
288 297
 	##use.ff=cnOptions[["use.ff"]]
289 298
 	outdir <- cnOptions[["outdir"]]
290
-	tmpdir <- cnOptions[["tmpdir"]]
291 299
 	if(missing(cdfName)) stop("cdfName is missing -- a valid cdfName is required.  See crlmm:::validCdfNames()")
292 300
 	platform <- whichPlatform(cdfName)
293 301
 	if(!(platform %in% c("affymetrix", "illumina"))){
... ...
@@ -454,18 +462,6 @@ whichPlatform <- function(cdfName){
454 462
 }
455 463
 
456 464
 
457
-##isValidCdfName <- function(cdfName, platform){
458
-##	chipList <- validCdfNames(platform)
459
-##	if(!(cdfName %in% chipList)){
460
-##		warning("cdfName must be one of the following: ",
461
-##			chipList)
462
-##	}
463
-##	result <- cdfName %in% chipList
464
-##	return(result)
465
-##}
466
-	
467
-	
468
-	
469 465
 # steps: quantile normalize hapmap: create 1m_reference_cn.rda object
470 466
 cnrma <- function(filenames, cdfName, sns, seed=1, verbose=FALSE, outdir){
471 467
 	if(missing(cdfName)) stop("must specify cdfName")
... ...
@@ -565,8 +561,7 @@ thresholdCopynumber <- function(object){
565 561
 ##
566 562
 ##}
567 563
 
568
-cnOptions <- function(tmpdir=tempdir(),
569
-		      outdir="./",
564
+cnOptions <- function(outdir="./",
570 565
 		      cdfName,
571 566
 		      crlmmFile="snpsetObject.rda",
572 567
 		      intensityFile="normalizedIntensities.rda",
... ...
@@ -581,14 +576,13 @@ cnOptions <- function(tmpdir=tempdir(),
581 576
 		      DF.PRIOR=50,
582 577
 		      bias.adj=FALSE,
583 578
 		      prior.prob=rep(1/4, 4),
584
-		      SNRmin=5,
579
+		      SNRmin=4,
585 580
 		      chromosome=1:24,
586 581
 		      seed=123,
587 582
 		      verbose=TRUE,
588 583
 		      GT.CONF.THR=0.99,
589 584
 		      PHI.THR=2^6,##used in nonpolymorphic fxn for training
590
-		      nAA.THR=5, ##used in nonpolymorphic fxn for training
591
-		      nBB.THR=5, ##used in nonpolymorphic fxn for training
585
+		      nHOM.THR=5, ##used in nonpolymorphic fxn for training
592 586
 		      MIN.NU=2^3,
593 587
 		      MIN.PHI=2^3,
594 588
 		      THR.NU.PHI=TRUE,
... ...
@@ -609,9 +603,8 @@ cnOptions <- function(tmpdir=tempdir(),
609 603
 ##		hmmOpts <- hmmOptions(...)
610 604
 ##	}
611 605
 	if(is.null(batch))
612
-		stop("batch must have the same length as the number of samples")
613
-	list(tmpdir=tmpdir,
614
-	     outdir=outdir,
606
+		stop("must specify batch -- should be the same length as the number of files")
607
+	list(outdir=outdir,
615 608
 	     cdfName=cdfName,
616 609
 	     crlmmFile=file.path(outdir, crlmmFile),
617 610
 	     intensityFile=file.path(outdir, intensityFile),
... ...
@@ -632,8 +625,7 @@ cnOptions <- function(tmpdir=tempdir(),
632 625
 	     seed=seed,
633 626
 	     verbose=verbose,
634 627
 	     PHI.THR=PHI.THR,
635
-	     nAA.THR=nAA.THR,
636
-	     nBB.THR=nBB.THR,
628
+	     nHOM.THR=nHOM.THR,
637 629
 	     MIN.NU=MIN.NU,
638 630
 	     MIN.PHI=MIN.PHI,
639 631
 	     THR.NU.PHI=THR.NU.PHI,
... ...
@@ -692,8 +684,8 @@ nonpolymorphic <- function(object, cnOptions, tmp.objects){
692 684
 		flags <- list(A=flagsA, B=flagsB)
693 685
 		return(flags)
694 686
 	}
695
-	nAA.THR <- cnOptions$nAA.THR
696
-	nBB.THR <- cnOptions$nBB.THR
687
+	nAA.THR <- cnOptions$nHOM.THR
688
+	nBB.THR <- cnOptions$nHOM.THR
697 689
 	PHI.THR <- cnOptions$PHI.THR
698 690
 	snpflags <- goodSnps(object, PHI.THR, tmp.objects, nAA.THR, nBB.THR)
699 691
 	flagsA <- snpflags$A
... ...
@@ -1115,21 +1107,14 @@ locationAndScale <- function(object, cnOptions, tmp.objects){
1115 1107
 	return(object)
1116 1108
 }
1117 1109
 
1118
-##coefs <- function(plateIndex, conf, MIN.OBS=3, envir, CONF.THR=0.99){
1119 1110
 coefs <- function(object, cnOptions, tmp.objects){
1120
-##	p <- plateIndex
1121 1111
 	batch <- unique(object$batch)
1122
-##	plates.completed <- envir[["plates.completed"]]
1123
-##	if(!plates.completed[p]) return()
1124
-##	CHR <- envir[["chrom"]]
1125 1112
 	CHR <- unique(chromosome(object))
1126
-##	plate <- envir[["plate"]]
1127 1113
 	muA <- tmp.objects[["muA"]]
1128 1114
 	muB <- tmp.objects[["muB"]]
1129 1115
 	vA <- tmp.objects[["vA"]]
1130 1116
 	vB <- tmp.objects[["vB"]]
1131 1117
 	Ns <- tmp.objects[["Ns"]]
1132
-##	uplate <- tmp.objects[["uplate"]]
1133 1118
 	if(CHR != 23){
1134 1119
 		IA <- muA[, 3:5]
1135 1120
 		IB <- muB[, 3:5]
... ...
@@ -1222,7 +1207,7 @@ polymorphic <- function(object, cnOptions, tmp.objects){
1222 1207
 	return(object)
1223 1208
 }
1224 1209
 
1225
-biasAdj <- function(object, cnOptions, tmp.objects){
1210
+posteriorProbability.snps <- function(object, cnOptions, tmp.objects=list()){
1226 1211
 	I <- isSnp(object)
1227 1212
 	gender <- object$gender
1228 1213
 	CHR <- unique(chromosome(object))
... ...
@@ -1300,6 +1285,18 @@ biasAdj <- function(object, cnOptions, tmp.objects){
1300 1285
 	posteriorProb[, , 2] <- hemDel
1301 1286
 	posteriorProb[, , 3] <- norm
1302 1287
 	posteriorProb[, , 4] <- amp
1288
+	return(list(tmp.objects, posteriorProb))
1289
+}
1290
+
1291
+biasAdj <- function(object, cnOptions, tmp.objects){
1292
+	gender <- object$gender
1293
+	CHR <- unique(chromosome(object))
1294
+	I <- isSnp(object)
1295
+	A <- A(object)[I, ]
1296
+	normal <- tmp.objects[["normal"]][I, ]
1297
+	results <- posteriorProbability.snps(object=object, cnOptions=cnOptions, tmp.objects=tmp.objects)
1298
+	posteriorProb <- results[[2]]
1299
+	tmp.objects <- results[[1]]
1303 1300
 	mostLikelyState <- apply(posteriorProb, c(1, 2), function(x) order(x, decreasing=TRUE)[1])
1304 1301
 	if(CHR == 23){
1305 1302
 		##so state index 3 is the most likely state for men and women
... ...
@@ -1430,35 +1427,6 @@ thresholdModelParams <- function(object, cnOptions){
1430 1427
 	return(object)
1431 1428
 }
1432 1429
 
1433
-
1434
-
1435
-##setMethod("update", "character", function(object, ...){
1436
-##	crlmmFile <- object
1437
-##	for(i in seq(along=crlmmFile)){
1438
-##		cat("Processing ", crlmmFile[i], "...\n")
1439
-##		load(crlmmFile[i])
1440
-##		crlmmSetList <- get("crlmmSetList")
1441
-##		if(length(crlmmSetList) == 3) next()  ##copy number object already present. 
1442
-##		if(!"chromosome" %in% fvarLabels(crlmmSetList[[1]])){
1443
-##			featureData(crlmmSetList[[1]]) <- addFeatureAnnotation(crlmmSetList)
1444
-##		} 
1445
-##		CHR <- unique(chromosome(crlmmSetList[[1]]))
1446
-##		if(length(CHR) > 1) stop("More than one chromosome in the object. This method requires one chromosome at a time.")		
1447
-##		if(CHR==24){
1448
-##			message("skipping chromosome 24")
1449
-##			next()
1450
-##		}
1451
-##		cat("----------------------------------------------------------------------------\n")
1452
-##		cat("-        Estimating copy number for chromosome", CHR, "\n")
1453
-##		cat("----------------------------------------------------------------------------\n")		
1454
-##		crlmmSetList <- update(crlmmSetList, CHR=CHR, ...)
1455
-##		save(crlmmSetList, file=crlmmFile[i])
1456
-##		rm(crlmmSetList); gc();
1457
-##	}
1458
-##})
1459
-
1460
-
1461
-
1462 1430
 computeCopynumber.SnpSuperSet <- function(object, cnOptions){
1463 1431
 ##	use.ff <- cnOptions[["use.ff"]]
1464 1432
 ##	if(!use.ff){
... ...
@@ -1478,16 +1446,6 @@ computeCopynumber.SnpSuperSet <- function(object, cnOptions){
1478 1446
 	return(object)
1479 1447
 }
1480 1448
 
1481
-## computeCopynumber.CNSet
1482
-##    tmp.objects <- instantiateObjects() ##vA, vB, (background mean), muA, muB (within genotype means), Ns
1483
-##    tmp.objects <- withinGenotypeMoments() ##compute vA, vB, muA, muB, Ns
1484
-##    object.batch <- locationAndScale()  ##calcualte tau2A, tau2B, sig2A, sig2B, corr for polymorphic probes
1485
-##    tmp.objects <- oneBatch()  ##impute muA, muB where necessary
1486
-##    object.batch <- coefs()    ##calls nuphiAllele
1487
-##            nuPhiAllele() ##updates nuA, nuA.se, phiA, phiA.se, nuB, ...
1488
-##    object.batch <- polymorphic()        assigns CA, CB
1489
-##    object.batch <- nonpolymorphic()     assigns CA
1490
-
1491 1449
 
1492 1450
 computeCopynumber.CNSet <- function(object, cnOptions){
1493 1451
 	CHR <- unique(chromosome(object))
... ...
@@ -1563,3 +1521,4 @@ computeCopynumber.CNSet <- function(object, cnOptions){
1563 1521
 
1564 1522
 
1565 1523
 
1524
+
... ...
@@ -166,3 +166,35 @@ ellipse.CNSet <- function(x, copynumber, batch, ...){
166 166
 
167 167
 
168 168
 
169
+setMethod("show", "CNSet", function(object){
170
+	callNextMethod(object)
171
+##	cat("emissionPr\n")
172
+##	cat("   array:", nrow(object), "features,", ncol(object), "samples,", dim(emissionPr(object))[3], "states\n")
173
+##	cat("   hidden states:\n")
174
+##	cat("      ", dimnames(emissionPr(object))[[3]], "\n")
175
+##	cat("   Missing values:", sum(is.na(emissionPr(object))), "\n")
176
+##	if(!all(is.na(emissionPr(object)))){
177
+##		cat("   minimum value:", min(emissionPr(object), na.rm=TRUE), "\n")
178
+##	} else  cat("   minimum value: NA (all missing)\n")
179
+##	cat("rangedData:  ")
180
+##	cat("    ", show(rangedData(object)), "\n")
181
+})
182
+##setMethod("rangedData", "CNSet", function(object) segmentData(object))
183
+##setReplaceMethod("rangedData", c("CNSet", "RangedData"), function(object, value){
184
+##	segmentData(object) <- value
185
+##})
186
+##
187
+##setMethod("segmentData", "CNSet", function(object) object@segmentData)
188
+##setReplaceMethod("segmentData", c("CNSet", "RangedData"), function(object, value){
189
+##	object@segmentData <- value
190
+##	object
191
+##})
192
+##
193
+##setMethod("emissionPr", "CNSet", function(object) object@emissionPr)
194
+##setReplaceMethod("emissionPr", c("CNSet", "array"), function(object, value){
195
+##	object@emissionPr <- value
196
+##	object
197
+##})
198
+##setMethod("start", "CNSet", function(x, ...) start(segmentData(x), ...))
199
+##setMethod("end", "CNSet", function(x, ...) end(segmentData(x), ...))
200
+##setMethod("width", "CNSet", function(x) width(segmentData(x)))
... ...
@@ -49,6 +49,7 @@ samples on SNP 5.0 made available via the \Rpackage{hapmapsnp5}
49 49
 package.
50 50
 
51 51
 <<crlmm>>=
52
+require(oligoClasses)
52 53
 library(crlmm)
53 54
 library(hapmapsnp5)
54 55
 path <- system.file("celFiles", package="hapmapsnp5")
... ...
@@ -43,211 +43,60 @@ options(prompt="R> ")
43 43
 CRLMM supports the following platforms:
44 44
 
45 45
 <<cdfname>>=
46
+library(oligoClasses)
46 47
 library(crlmm)
47 48
 crlmm:::validCdfNames()
48 49
 @ 
49 50
 
50
-<<Rsge, eval=TRUE, echo=FALSE>>=
51
-##If multiple cluster nodes are available
52
-library(Rsge)
53
-@ 
54
-
55 51
 \paragraph{Preprocess and genotype.}
56 52
 
57
-Specify the coordinates of Affymetrix cel files and where to put
58
-intermediate files generated during the course of preprocessing / copy
59
-number estimation.
53
+Provide the complete path for the filenames:
60 54
 
61 55
 <<celfiles>>=
62
-myPath <- "/thumper/ctsa/snpmicroarray/hapmap/raw/affy/1m"
63
-celFiles <- list.celfiles(myPath, full.names=TRUE, pattern=".CEL")
64
-outdir <- "/thumper/ctsa/snpmicroarray/rs/data/hapmap/1m/affy"
65
-#outdir <- "/thumper/ctsa/snpmicroarray/rs/data/hapmap/1m/tmp"
66
-if(!file.exists(outdir)) dir.create(outdir)
67
-@ 
68
-
69
-Split the list of files by plate.  For this HapMap data, the different
70
-ancestries were run on separate plates.
56
+celFiles <- list.celfiles("/thumper/ctsa/snpmicroarray/hapmap/raw/affy/1m", full.names=TRUE, pattern=".CEL")
57
+@ 
58
+
59
+While genotyping with crlmm can be performed using a small number of
60
+samples, copy number estimation requires at least 10 samples --
61
+preferably all of the samples that were processed together on the same
62
+plate.  Crlmm does not use a reference dataset when estimating model
63
+parameters because of large batch effects \citep{Scharpf2009}. The
64
+quantile normalization performed as part of the preprocessing of the raw
65
+data is insufficient for removing batch effects.  Processing a reference
66
+dataset, such as HapMap samples, along with the experimental data will
67
+not improve copy number estimation for the experimental dataset, and
68
+should not be used as a means to increase the sample size. Furthermore,
69
+processing a reference dataset without acknowledging that these samples
70
+were derived from a different batch can result in incorrect copy number
71
+estimates in both the experimental and reference datasets.  The
72
+appropriate way to acknowledge batch is to supply the batch name for
73
+each sample to be processed in the argument to \Robject{cnOptions}:
71 74
 
72
-<<plate>>=
73
-plate <- substr(basename(celFiles), 13, 13)
74
-table(plate)
75
-celFiles <- split(celFiles, plate)
75
+<<celfiles>>=
76
+cnOpts <- cnOptions(cdfName="genomewidesnp6",
77
+		    outdir="/thumper/ctsa/snpmicroarray/hapmap/raw/affy/1m",
78
+		    batch=substr(basename(celFiles), 13, 13))
79
+str(cnOpts)
80
+stopifnot(length(cnOpts$batch) == length(celFiles))
81
+names(cnOpts$batch) <- basename(celFiles)
76 82
 @ 
77 83
 
78
-
79 84
 \noindent The next code chunk quantile normalizes the samples to a
80
-target reference distribution and uses the crlmm algorithm to genotype.
81
-See \cite{Carvalho2007a} for details regarding the crlmm genotyping
82
-algorithm.
85
+target reference distribution, uses the crlmm algorithm to genotype, and
86
+then estimates the copy number for each batch. Currently processing the
87
+270 HapMap cel files will require approximately 20G of RAM. We are
88
+working on methods to reduce the memory footprint.
83 89
 
84 90
 <<preprocessAndGenotype, eval=TRUE, echo=TRUE>>=
85
-for(i in seq(along=celFiles)){
86
-	platedir <- file.path(outdir, names(celFiles)[i])
87
-	if(!file.exists(platedir)) dir.create(platedir)
88
-	crlmmWrapper(celFiles[[i]],##[1:15], 
89
-		     cdfName="genomewidesnp6",
90
-		     save.it=FALSE, 
91
-		     load.it=TRUE,
92
-		     intensityFile=file.path(platedir, "normalizedIntensities.rda"),
93
-		     crlmmFile=file.path(platedir, "snpsetObject.rda"))
94
-}
95
-q("no")
96
-@ 
97
-
98
-<<preprocessAndGenotype, eval=FALSE, echo=FALSE>>=
99
-for(i in seq(along=celFiles)){
100
-	platedir <- file.path(outdir, names(celFiles)[i])
101
-	if(!file.exists(platedir)) dir.create(platedir)
102
-	sge.submit(crlmmWrapper,
103
-		   filenames=celFiles[[i]],##[1:15], 
104
-		   cdfName="genomewidesnp6",
105
-		   save.it=FALSE, 
106
-		   load.it=TRUE,
107
-		   intensityFile=file.path(platedir, "normalizedIntensities.rda"),
108
-		   crlmmFile=file.path(platedir, "snpsetObject.rda"),
109
-		   global.savelist=c("platedir"),
110
-		   packages="crlmm")
111
-}
112
-##platedirs <- sapply(names(celFiles), function(x) file.path(outdir, x))
113
-##intensityFiles <- file.path(platedirs, "normalizedIntensities.rda")
114
-####crlmmFiles <- file.path(platedirs, "snpsetObject.rda")
115
-####myF <- function(celFiles, intensityFile, crlmmFile, ...){
116
-####	crlmmWrapper(celFiles, intensityFile=intensityFile,
117
-####		     crlmmFile=crlmmFile,
118
-##sge.apply(celFiles, 
119
-##	  crlmmWrapper,
120
-##	  cdfName="genomewidesnp6",
121
-##	  save.it=FALSE, 
122
-##	  load.it=TRUE,
123
-##	  intensityFile=,
124
-##	  crlmmFile=file.path(platedir, "snpsetObject.rda"))
125
-q("no")
91
+crlmmCopynumber(celFiles, cnOpts)
126 92
 @ 
127 93
 
128
-As a result of the above wrapper, the following R objects are now
129
-created for each plate:
94
+The following R objects are created from crlmmCopynumber:
130 95
 
131 96
 <<crlmmSetListObjects>>=
132
-fns <- list.files(platedir, pattern="crlmmSetList", full.name=TRUE)
97
+fns <- list.files(cnOpts[["outdir"]], pattern="cnSet", full.name=TRUE)
133 98
 @ 
134 99
 
135
-To estimate copy number, simply run the update function on this vector
136
-of filenames.  If multiple nodes are available for computation, one
137
-could easily split this job.
138
-
139
-<<copynumber, echo=FALSE, eval=TRUE>>=
140
-for(i in seq(along=celFiles)){
141
-	fns <- list.files(file.path(outdir, names(celFiles[i])), 
142
-			  pattern="crlmmSetList", full.names=TRUE)
143
-	update(fns)
144
-}
145
-@
146
-
147
-<<copynumber.sge, echo=FALSE, eval=TRUE>>=
148
-for(i in seq(along=celFiles)){
149
-	fns <- list.files(file.path(outdir, names(celFiles[i])), 
150
-			  pattern="crlmmSetList", full.names=TRUE)
151
-	sge.parLapply(X=fns, FUN=update, packages="crlmm", cluster=FALSE)
152
-}
153
-@
154
-
155
-\paragraph{Locus- and allele-specific estimates of copy number.}
156
-
157
-Load the object for chromosome 22 and compute copy number:
158
-
159
-<<crlmmList-show>>=
160
-CHR <- 22
161
-if(!exists("crlmmSetList")) load(file.path(outdir, paste("crlmmSetList_", CHR, ".rda", sep="")))
162
-show(crlmmSetList)
163
-crlmmSetList <- update(crlmmSetList)
164
-show(crlmmSetList)
165
-@ 
166
-
167
-See the help file for \Robject{computeCopynumber} for arguments to
168
-\Robject{update}.  Provided that the assumption of integer copy number
169
-is reasonable, one can fit a hidden Markov model to the locus-level
170
-estimates of copy number and uncertainty.
171
-
172
-\paragraph{A hidden Markov model.}
173
-
174
-Emission probabilities for the hidden markov model can be computed from
175
-the copy number prediction regions based on bivariate normal scatter
176
-plots of the log A versus log B intensities.  (By contrast, a
177
-nonparamemtric segmentation would be applied to intensities on the scale
178
-of copy number.)
179
-
180
-<<emissionProbabilities>>=
181
-library(VanillaICE)
182
-copyNumberStates <- 0:5
183
-if(!exists("emission.cn")){
184
-	emission.cn <- suppressWarnings(crlmm:::computeEmission(crlmmSetList, copyNumberStates))
185
-	dim(emission.cn)
186
-}
187
-@ 
188
-
189
-*Before smoothing the estimates as a function of physical position by
190
-segmentation hidden Markov models, one should consider how to handle
191
-outliers.  In particular, samples with low signal to noise ratios are
192
-likely to have many copy number outliers. See suggested visualizations.
193
-
194
-Initial state probabilities and transition probabilities for the HMM:
195
-
196
-<<probs>>=
197
-initialP <- rep(1/length(copyNumberStates), length(copyNumberStates))
198
-tau <- transitionProbability(chromosome=chromosome(crlmmSetList),
199
-			     position=position(crlmmSetList),
200
-			     TAUP=1e8)
201
-@ 
202
-
203
-The viterbi algorithm is used to identify the sequence of states that
204
-maximizes the likelihood:
205
-
206
-<<hmm>>=
207
-if(!exists("hmmPredictions")){
208
-hmmPredictions <- viterbi(emission=emission.cn,
209
-			  initialStateProbs=log(initialP), 
210
-			  tau=tau[, "transitionPr"],
211
-			  arm=tau[, "arm"], 
212
-			  normalIndex=3,
213
-			  normal2altered=0.01,
214
-			  altered2normal=1,
215
-			  altered2altered=0.001)	    
216
-}
217
-table(as.integer(hmmPredictions)-1)
218
-brks <- breaks(x=hmmPredictions, states=copyNumberStates, 
219
-	       position=tau[, "position"],
220
-	       chromosome=tau[, "chromosome"])
221
-str(brks)
222
-@ 
223
-
224
-\paragraph{Circular binary segmentation}
225
-
226
-<<setup>>=
227
-library(DNAcopy)
228
-library(genefilter)
229
-ratioset <- crlmmSetList[[3]]
230
-meds <- rowMedians(copyNumber(ratioset), na.rm=TRUE)
231
-ratios <- log2(copyNumber(ratioset)) - log2(meds)
232
-ratioset <- new("SnpCopyNumberSet",
233
-		copyNumber=ratios,
234
-		phenoData=phenoData(ratioset),
235
-		featureData=featureData(ratioset),
236
-		annotation=annotation(ratioset))
237
-@ 
238
-
239
-The following code chunk (not evaluated) takes a while to run:
240
-
241
-<<cbs, eval=FALSE>>=
242
-CNA.object <- CNA(copyNumber(ratioset), 
243
-		  chromosome(ratioset),
244
-		  position(ratioset),
245
-		  data.type="logratio",
246
-		  sampleid=sampleNames(ratioset))
247
-smoothed.CNA.object <- smooth.CNA(CNA.object)
248
-segment.cna <- segment(smoothed.CNA.object, verbose=1)
249
-save(segment.cna, file=file.path(outdir, paste("segment.cna_", CHR, ".rda", sep="")))
250
-@ 
251 100
 
252 101
 \section{Accessors}
253 102
 
... ...
@@ -258,69 +107,59 @@ save(segment.cna, file=file.path(outdir, paste("segment.cna_", CHR, ".rda", sep=
258 107
 An object of class \Robject{ABset} is stored in the first element of the
259 108
 \Robject{crlmmSetList} object. The following accessors may be of use:
260 109
 
261
-Accessors for the quantile normalized intensities for the A allele:
110
+Accessors for the quantile normalized intensities for the A allele at
111
+polymorphic loci:
262 112
 
263 113
 <<accessors>>=
264
-a <- A(crlmmSetList)
114
+a <- A(cnSet)[isSnp(cnSet), ]
265 115
 dim(a)
266 116
 @ 
267 117
 
268
-The quantile normalized intensities for the nonpolymorphic probes are
269
-also stored in the 'A' assay data element.  To retrieve the quantile
270
-normalized intensities for the A allele only at polymorphic loci:
271
-
272
-<<A.polymorphic>>=
273
-a.snps <- A(crlmmSetList[snpIndex(crlmmSetList), ])
274
-dim(a.snps)
275
-@ 
118
+The quantile-normalized intensities for nonpolymorphic loci are obtained by:
276 119
 
277
-\noindent For the nonpolymorphic loci:
278
-<<A.nonpolymorphic>>=
279
-a.nps <- A(crlmmSetList[cnIndex(crlmmSetList), ])
280
-dim(a.nps)
120
+<<>>=
121
+npIntensities <- A(cnSet)[!isSnp(cnSet), ]
281 122
 @ 
282 123
 
283 124
 Quantile normalized intensities for the B allele at polymorphic loci:
284 125
 
285 126
 <<B.polymorphic>>=
286
-b.snps <- B(crlmmSetList[snpIndex(crlmmSetList), ])
127
+b.snps <- B(cnSet[isSnp(cnSet), ])
287 128
 @ 
288 129
 
289 130
 Note that NAs are recorded in the 'B' assay data element for
290 131
 nonpolymorphic loci:
291 132
 
292 133
 <<B.NAs>>=
293
-all(is.na(B(crlmmSetList[cnIndex(crlmmSetList), ])))
134
+all(is.na(B(cnSet[!isSnp(cnSet), ])))
294 135
 @ 
295 136
 
296 137
 \paragraph{\Robject{SnpSet}: Genotype calls and confidence scores}
297 138
 
298 139
 Genotype calls:
299 140
 <<genotypes>>=
300
-genotypes <- calls(crlmmSetList)
141
+genotypes <- snpCall(cnSet)
301 142
 @ 
302 143
 Confidence scores of the genotype calls:
303 144
 <<confidenceScores>>=
304
-genotypeConf <- confs(crlmmSetList)
145
+genotypeConf <- confs(cnSet)
305 146
 @ 
306 147
 
307 148
 \paragraph{\Robject{CopyNumberSet}: allele-specific copy number}
308 149
 
309 150
 Allele-specific copy number at polymorphic loci:
310 151
 <<ca>>=
311
-ca <- CA(crlmmSetList[snpIndex(crlmmSetList), ])
152
+ca <- CA(cnSet[isSnp(cnSet), ])
312 153
 @ 
313 154
 
314 155
 Total copy number at nonpolymorphic loci:
315 156
 <<ca>>=
316
-cn.nonpolymorphic <- CA(crlmmSetList[cnIndex(crlmmSetList), ])
317
-range(cn.nonpolymorphic, na.rm=TRUE)
318
-median(cn.nonpolymorphic, na.rm=TRUE)
157
+cn.nonpolymorphic <- CA(cnSet[!isSnp(cnSet), ])
319 158
 @ 
320 159
 
321 160
 Total copy number at both polymorphic and nonpolymorphic loci:
322 161
 <<totalCopynumber>>=
323
-cn <- copyNumber(crlmmSetList)
162
+cn <- copyNumber(cnSet)
324 163
 @ 
325 164
 
326 165
 \subsection{Other accessors}
... ...
@@ -347,7 +186,7 @@ fvarLabels(crlmmSetList[[3]])
347 186
 A histogram of the signal to noise ratio for the HapMap samples:
348 187
 
349 188
 <<plotSnr, fig=TRUE, include=FALSE>>=
350
-hist(crlmmSetList[[2]]$SNR, xlab="SNR", main="")
189
+hist(cnSet$SNR, xlab="SNR", main="")
351 190
 @ 
352 191
 
353 192
 \begin{figure}
... ...
@@ -367,12 +206,12 @@ one would have 70+ files in a given batch. Here we make a table of date
367 206
 versus ancestry:
368 207
 
369 208
 <<specifyBatch, eval=FALSE, echo=FALSE>>=
370
-sns <- sampleNames(crlmmResults)
209
+sns <- sampleNames(cnSet)
371 210
 sns[1]
372 211
 plate <- substr(basename(sns), 13, 13)
373 212
 table(plate)
374
-table(format(as.POSIXlt(protocolData(crlmmResults)[["ScanDate"]]), "%d %b %Y"), plate)
375
-dts <- protocolData(crlmmResults)[["ScanDate"]]
213
+table(format(as.POSIXlt(protocolData(cnSet)[["ScanDate"]]), "%d %b %Y"), plate)
214
+dts <- protocolData(cnSet)[["ScanDate"]]
376 215
 @ 
377 216
 
378 217
 As all of these samples were run on the first week of March, we would
... ...
@@ -401,69 +240,70 @@ additional robustness to this assumption.  Set the \Robject{bias.adj}
401 240
 argument to \Robject{TRUE}:
402 241
 
403 242
 <<biasAdjustment, eval=FALSE>>=
404
-update(crlmmSetList, CHR=22, bias.adj=TRUE)
405
-@ 
406
-
407
-The following code chunk calculates the frequency of amplifications and
408
-deletions at each locus. Shaded regions above the zero line in Figure
409
-\ref{fig:hmm_hapmap} display the frequency of amplifications, whereas
410
-shaded regions below the zero line graphically display the frequency of
411
-hemizygous or homozygous deletions.
412
-
413
-<<hmm_hapmap, fig=TRUE, include=FALSE, width=8, height=7>>=
414
-require(SNPchip)
415
-library(RColorBrewer)
416
-numberUp <- rowSums(hmmPredictions > 3, na.rm=TRUE)
417
-numberDown <- -rowSums(hmmPredictions < 3, na.rm=TRUE)
418
-poly.cols <- brewer.pal(7, "Accent")
419
-alt.brks <- brks[brks[, "state"] != "copy.number_2", ]
420
-op <- par(ask=FALSE)
421
-ylim <- c(min(numberDown)-5, max(numberUp)+5)
422
-xlim <- c(10*1e6, max(position(crlmmSetList)))
423
-plot(position(crlmmSetList), rep(0, nrow(crlmmSetList[[1]])),
424
-     type="n", xlab="Physical position (Mb)",
425
-     ylim=ylim,
426
-     xlim=xlim,
427
-     ylab="frequency", main="Chr 22",
428
-     xaxt="n",
429
-     xaxs="r")
430
-axis(1, at=pretty(xlim), labels=pretty(xlim)/1e6)
431
-polygon(x=c(position(crlmmSetList), rev(position(crlmmSetList))),
432
-	y=c(rep(0, nrow(crlmmSetList[[1]])), rev(numberUp)),
433
-	col=poly.cols[3], border=poly.cols[3])
434
-polygon(x=c(position(crlmmSetList), rev(position(crlmmSetList))),
435
-	y=c(rep(0, nrow(crlmmSetList[[1]])), rev(numberDown)),
436
-	col=poly.cols[5], border=poly.cols[5])
437
-##plotCytoband(22, xlim=xlim, new=FALSE,
438
-##	     label.cytoband=FALSE,
439
-##	     cytoband.ycoords=c(-10, -8), xaxs="r")
440
-
441
-medLength <- round(median(alt.brks[, "nbases"]), 2)
442
-medMarkers <- median(alt.brks[, "nprobes"])
443
-sdMarkers <- round(mad(alt.brks[, "nprobes"]), 2)
444
-sdsLength <- round(mad(alt.brks[, "nbases"]), 2)
445
-legend("topright",
446
-       bty="n",
447
-       legend=c(paste("median length:", medLength, "(bp)"),
448
-       paste("MAD length:", sdsLength, "(bp)"),
449
-       paste("median # markers:", medMarkers),
450
-       paste("MAD # markers:", sdMarkers)),
451
-       cex=0.8, ncol=2)
452
-legend("topleft",
453
-       fill=poly.cols[c(3, 5)],
454
-       legend=c("amplifications", "deletions"), bty="n")
455
-par(op)
456
-gc()
243
+cnOpts[["bias.adj"]] <- TRUE
244
+crlmmCopynumber(celFiles, cnOpts)
457 245
 @ 
458
-
459
-\begin{figure}
460
-  \centering
461
-  \includegraphics[width=\textwidth]{copynumber-hmm_hapmap}
462
-  \caption{\label{fig:hmm_hapmap} The frequency of amplifications in the
463
-    hapmap samples is displayed above the zero line.  The frequency of
464
-    hemizygous or homozygous deletions are displayed below the zero
465
-    line.}
466
-\end{figure}
246
+%
247
+%The following code chunk calculates the frequency of amplifications and
248
+%deletions at each locus. Shaded regions above the zero line in Figure
249
+%\ref{fig:hmm_hapmap} display the frequency of amplifications, whereas
250
+%shaded regions below the zero line graphically display the frequency of
251
+%hemizygous or homozygous deletions.
252
+%
253
+%<<hmm_hapmap, fig=TRUE, include=FALSE, width=8, height=7>>=
254
+%require(SNPchip)
255
+%library(RColorBrewer)
256
+%numberUp <- rowSums(hmmPredictions > 3, na.rm=TRUE)
257
+%numberDown <- -rowSums(hmmPredictions < 3, na.rm=TRUE)
258
+%poly.cols <- brewer.pal(7, "Accent")
259
+%alt.brks <- brks[brks[, "state"] != "copy.number_2", ]
260
+%op <- par(ask=FALSE)
261
+%ylim <- c(min(numberDown)-5, max(numberUp)+5)
262
+%xlim <- c(10*1e6, max(position(crlmmSetList)))
263
+%plot(position(crlmmSetList), rep(0, nrow(crlmmSetList[[1]])),
264
+%     type="n", xlab="Physical position (Mb)",
265
+%     ylim=ylim,
266
+%     xlim=xlim,
267
+%     ylab="frequency", main="Chr 22",
268
+%     xaxt="n",
269
+%     xaxs="r")
270
+%axis(1, at=pretty(xlim), labels=pretty(xlim)/1e6)
271
+%polygon(x=c(position(crlmmSetList), rev(position(crlmmSetList))),
272
+%	y=c(rep(0, nrow(crlmmSetList[[1]])), rev(numberUp)),
273
+%	col=poly.cols[3], border=poly.cols[3])
274
+%polygon(x=c(position(crlmmSetList), rev(position(crlmmSetList))),
275
+%	y=c(rep(0, nrow(crlmmSetList[[1]])), rev(numberDown)),
276
+%	col=poly.cols[5], border=poly.cols[5])
277
+%##plotCytoband(22, xlim=xlim, new=FALSE,
278
+%##	     label.cytoband=FALSE,
279
+%##	     cytoband.ycoords=c(-10, -8), xaxs="r")
280
+%
281
+%medLength <- round(median(alt.brks[, "nbases"]), 2)
282
+%medMarkers <- median(alt.brks[, "nprobes"])
283
+%sdMarkers <- round(mad(alt.brks[, "nprobes"]), 2)
284
+%sdsLength <- round(mad(alt.brks[, "nbases"]), 2)
285
+%legend("topright",
286
+%       bty="n",
287
+%       legend=c(paste("median length:", medLength, "(bp)"),
288
+%       paste("MAD length:", sdsLength, "(bp)"),
289
+%       paste("median # markers:", medMarkers),
290
+%       paste("MAD # markers:", sdMarkers)),
291
+%       cex=0.8, ncol=2)
292
+%legend("topleft",
293
+%       fill=poly.cols[c(3, 5)],
294
+%       legend=c("amplifications", "deletions"), bty="n")
295
+%par(op)
296
+%gc()
297
+%@ 
298
+%
299
+%\begin{figure}
300
+%  \centering
301
+%  \includegraphics[width=\textwidth]{copynumber-hmm_hapmap}
302
+%  \caption{\label{fig:hmm_hapmap} The frequency of amplifications in the
303
+%    hapmap samples is displayed above the zero line.  The frequency of
304
+%    hemizygous or homozygous deletions are displayed below the zero
305
+%    line.}
306
+%\end{figure}
467 307
 
468 308
 \clearpage
469 309
 \paragraph{One sample at a time: locus-level estimates}
... ...
@@ -473,12 +313,12 @@ versus copy number (vertical axis) for the first sample.
473 313
 
474 314
 <<oneSample, fig=TRUE, width=8, height=4, include=FALSE>>=
475 315
 par(las=1)
476
-plot(position(crlmmSetList), copyNumber(crlmmSetList)[, 1], pch=".", cex=2, xaxt="n", col="grey20", ylim=c(0,6), 
316
+plot(position(cnSet), copyNumber(cnSet)[, 1], pch=".", cex=2, xaxt="n", col="grey20", ylim=c(0,6), 
477 317
      ylab="copy number", xlab="physical position (Mb)",
478
-     main=paste(sampleNames(crlmmSetList)[1], ", CHR:", unique(chromosome(crlmmSetList))))
479
-points(position(crlmmSetList)[cnIndex(crlmmSetList)], copyNumber(crlmmSetList)[cnIndex(crlmmSetList), 1],
318
+     main=paste(sampleNames(cnSet)[1], ", CHR:", unique(chromosome(cnSet))))
319
+points(position(cnSet)[cnIndex(cnSet)], copyNumber(cnSet)[cnIndex(cnSet), 1],
480 320
        pch=".", cex=2, col="lightblue")
481
-axis(1, at=pretty(range(position(crlmmSetList))), labels=pretty(range(position(crlmmSetList)))/1e6)
321
+axis(1, at=pretty(range(position(cnSet))), labels=pretty(range(position(cnSet)))/1e6)
482 322
 @ 
483 323
 
484 324
 <<idiogram, eval=FALSE, echo=FALSE>>=
... ...
@@ -501,20 +341,20 @@ estimates and overlays the predictions from the hidden markov model.
501 341
 ask <- FALSE
502 342
 op <- par(mfrow=c(3, 1), las=1, mar=c(1, 4, 1, 1), oma=c(3, 1, 1, 1), ask=ask)
503 343
 ##Put fit on the copy number scale
504
-cns <- copyNumber(crlmmSetList)
344
+cns <- copyNumber(cnSet)
505 345
 cnState <- hmmPredictions - as.integer(1)
506
-xlim <- c(10*1e6, max(position(crlmmSetList)))
346
+xlim <- c(10*1e6, max(position(cnSet)))
507 347
 cols <- brewer.pal(8, "Dark2")[1:4]
508 348
 for(j in 1:3){
509
-	plot(position(crlmmSetList), cnState[, j], pch=".", col=cols[2], xaxt="n", 
349
+	plot(position(cnSet), cnState[, j], pch=".", col=cols[2], xaxt="n", 
510 350
 	     ylab="copy number", xlab="Physical position (Mb)", type="s", lwd=2,
511 351
 	     ylim=c(0,6), xlim=xlim)
512
-	points(position(crlmmSetList), cns[, j], pch=".", col=cols[3])
513
-	lines(position(crlmmSetList), cnState[,j], lwd=2, col=cols[2])
514
-	axis(1, at=pretty(position(crlmmSetList)), 
515
-	     labels=pretty(position(crlmmSetList))/1e6)
352
+	points(position(cnSet), cns[, j], pch=".", col=cols[3])
353
+	lines(position(cnSet), cnState[,j], lwd=2, col=cols[2])
354
+	axis(1, at=pretty(position(cnSet)), 
355
+	     labels=pretty(position(cnSet))/1e6)
516 356
 	abline(h=c(1,3), lty=2, col=cols[1])	
517
-	legend("topright", bty="n", legend=sampleNames(crlmmSetList)[j])
357
+	legend("topright", bty="n", legend=sampleNames(cnSet)[j])
518 358
 	legend("topleft", lty=1, col=cols[2], legend="copy number state",
519 359
 	       bty="n", lwd=2)
520 360
 	plotCytoband(CHR, cytoband.ycoords=c(5, 5.2), new=FALSE,
... ...
@@ -545,12 +385,12 @@ colors <- c("red", "blue", "green3")
545 385
 cex <- 0.6
546 386
 par(mfrow=c(3,3), las=1, pty="s", ask=FALSE, mar=c(2, 2, 2, 2), oma=c(2, 2, 1, 1))
547 387
 ##plot 9 at a time
548
-indices <- split(snpIndex(crlmmSetList), rep(1:length(snpIndex(crlmmSetList)), each=9, length.out=length(snpIndex(crlmmSetList))))
388
+indices <- split(snpIndex(cnSet), rep(1:length(snpIndex(cnSet)), each=9, length.out=length(snpIndex(cnSet))))
549 389
 ##par(ask=FALSE)
550 390
 j <- 1
551 391
 for(i in indices[[j]]){
552
-	gt <- calls(crlmmSetList)[i, ]
553
-	plot(crlmmSetList[i, ], 
392
+	gt <- calls(cnSet)[i, ]
393
+	plot(cnSet[i, ], 
554 394
 	     pch=pch, 
555 395
 	     col=colors[gt], 
556 396
 	     bg=colors[gt], cex=cex,
... ...
@@ -577,7 +417,7 @@ for(i in indices[[j]]){
577 417
 require(RColorBrewer)
578 418
 library(ellipse)
579 419
 greens <- brewer.pal(9, "Greens")
580
-J <- split(1:ncol(crlmmSetList), batch(crlmmSetList))
420
+J <- split(1:ncol(cnSet), batch(cnSet))
581 421
 colors <- c("red", "blue", "green3")
582 422
 cex <- 0.6
583 423
 colors <- c("blue", greens[8], "red")
... ...
@@ -588,33 +428,33 @@ lwd <- 2
588 428
 ##pdf("figures/snp22plots%02d.pdf", width=600, height=600)
589 429
 ask <- FALSE
590 430
 par(mfrow=c(3,3), las=1, pty="s", ask=ask, mar=c(2, 2, 2, 2), oma=c(2, 2, 1, 1))
591
-indices <- split(snpIndex(crlmmSetList), rep(1:length(snpIndex(crlmmSetList)), each=9, length.out=length(snpIndex(crlmmSetList))))
431
+indices <- split(snpIndex(cnSet), rep(1:length(snpIndex(cnSet)), each=9, length.out=length(snpIndex(cnSet))))
592 432
 ##for(j in seq(along=indices)[1:10]){
593 433
 j <- 1
594 434
 	cat(j, "\n")
595 435
 	k <- 1
596 436
 	for(i in indices[[j]]){
597
-		gt <- calls(crlmmSetList)[i, ]
437
+		gt <- calls(cnSet)[i, ]
598 438
 		pch <- as.character(gt)
599 439
 		cex <- 0.9
600
-		plot(crlmmSetList[i, ], 
440
+		plot(cnSet[i, ], 
601 441
 		     pch=pch, 
602 442
 		     col=pch.col[gt],
603 443
 		     cex=cex,
604 444
 		     xlim=xlim, ylim=ylim,
605 445
 		     type="n")
606 446
 		if(plotpoints){
607
-			for(b in seq(along=unique(batch(crlmmSetList)))){
608
-				points(crlmmSetList[i, J[[b]]], 
447
+			for(b in seq(along=unique(batch(cnSet)))){
448
+				points(cnSet[i, J[[b]]], 
609 449
 				       pch=pch, 
610 450
 				       col=colors[b], bg=colors[b], cex=cex,
611 451
 				       xlim=xlim, ylim=ylim)
612 452
 			}
613 453
 		}
614
-		for(b in seq(along=unique(batch(crlmmSetList)))){
615
-			ellipse(crlmmSetList[i, J[[b]]], copynumber=2, col=colors[b], lwd=lwd)
454
+		for(b in seq(along=unique(batch(cnSet)))){
455
+			ellipse(cnSet[i, J[[b]]], copynumber=2, col=colors[b], lwd=lwd)
616 456
 		}
617
-		##legend("bottomright", bty="n", legend=featureNames(crlmmSetList)[i])
457
+		##legend("bottomright", bty="n", legend=featureNames(cnSet)[i])
618 458
 		if(k == 1) {
619 459
 			legend("bottomleft", bty="n", fill=colors, legend=c("CEPH", "Yoruba", "Asian"))	
620 460
 			mtext("A", 1, outer=TRUE, line=1)
... ...
@@ -645,11 +485,11 @@ adjusted by passing the argument \Robject{CONF.THR} to the
645 485
 
646 486
 <<linearModel, fig=TRUE,include=FALSE, eval=FALSE, echo=FALSE>>=
647 487
 cols <- brewer.pal(7, "Accent")[c(1, 4, 7)]
648
-i <- match("SNP_A-8608180", featureNames(crlmmSetList))
649
-B <- batch(crlmmSetList) == unique(batch(crlmmSetList))[1]
650
-gt.conf <- confs(crlmmSetList[i, B])
488
+i <- match("SNP_A-8608180", featureNames(cnSet))
489
+B <- batch(cnSet) == unique(batch(cnSet))[1]
490
+gt.conf <- confs(cnSet[i, B])
651 491
 op <- par(mfrow=c(2, 1), las=1, ask=FALSE, mar=c(3.5, 4, 0.5, 0.5), oma=c(0, 0, 2, 0))
652
-crlmm:::boxplot(crlmmSetList[i, B][[1]])
492
+crlmm:::boxplot(cnSet[i, B][[1]])
653 493
 par(op)
654 494
 @ 
655 495
 
... ...
@@ -45,6 +45,7 @@ are assumed to be in the current directory. First, we identify the
45 45
 files and run \Rmethod{crlmm}. The results will be saved to the
46 46
 variable \Robject{crlmmResult}.
47 47
 <<lkd>>=
48
+library(oligoClasses}
48 49
 library(crlmm)
49 50
 celFiles <- list.celfiles()
50 51
 celFiles[1:4]
... ...
@@ -8,8 +8,6 @@
8 8
 \alias{B,AlleleSet-method}
9 9
 \alias{B<-}
10 10
 \alias{B<-,AlleleSet,matrix-method}
11
-\alias{isSnp}
12
-\alias{isSnp,AlleleSet-method}
13 11
 \title{Indicator for polymorphic probes}
14 12
 \description{
15 13
   
... ...
@@ -23,12 +21,11 @@
23 21
 \usage{
24 22
 A(object)
25 23
 B(object)
26
-isSnp(object)
27 24
 }
28 25
 \arguments{
29 26
   \item{object}{AlleleSet object}
30 27
 }
31 28
 \value{
32
-  an indicator: 0=nonpolymorphic, 1=polymorphic
29
+	matrix of normalized intensities
33 30
 }
34 31
 \keyword{manip}
35 32
deleted file mode 100644
... ...
@@ -1,69 +0,0 @@
1
-\name{CNSet-class}
2
-\Rdversion{1.1}
3
-\docType{class}
4
-\alias{CNSet-class}
5
-\alias{CA}
6
-\alias{CA<-}
7
-\alias{CB}
8
-\alias{CB<-}
9
-\alias{CA,CNSet-method}
10
-\alias{CA<-,CNSet,matrix-method}
11
-\alias{CB,CNSet-method}
12
-\alias{CB<-,CNSet,matrix-method}
13
-\alias{computeHmm}
14
-\alias{computeHmm,CNSet-method}
15
-\alias{copyNumber,CNSet-method}
16
-\alias{ellipse,CNSet-method}
17
-\alias{emissionPr}
18
-\alias{emissionPr<-}
19
-\alias{emissionPr,CNSet-method}
20
-\alias{emissionPr<-,CNSet,array-method}
21
-\alias{segmentData}
22
-\alias{segmentData<-}
23
-\alias{segmentData,CNSet-method}
24
-\alias{segmentData<-,CNSet,RangedData-method}
25
-
26
-\title{Container for allele-specific copy number and lower level SNP /
27
-  NP summaries}
28
-\description{Container for allele-specific copy number and lower level SNP /
29
-  NP summaries}
30
-\section{Objects from the Class}{
31
-Objects can be created by calls of the form \code{new("CNSet", phenoData, featureData, annotation, experimentData, protocolData, calls, callsConfidence, senseThetaA, senseThetaB, CA, CB, ...)}.
32
-}
33
-\section{Slots}{
34
-  \describe{
35
-    \item{\code{assayData}:}{Object of class \code{"AssayData"} ~~ }
36
-    \item{\code{phenoData}:}{Object of class \code{"AnnotatedDataFrame"} ~~ }
37
-    \item{\code{featureData}:}{Object of class \code{"AnnotatedDataFrame"} ~~ }
38
-    \item{\code{experimentData}:}{Object of class \code{"MIAME"} ~~ }
39
-    \item{\code{annotation}:}{Object of class \code{"character"} ~~ }
40
-    \item{\code{protocolData}:}{Object of class \code{"AnnotatedDataFrame"} ~~ }
41
-    \item{\code{.__classVersion__}:}{Object of class \code{"Versions"} ~~ }
42
-  }
43
-}
44
-\section{Extends}{
45
-Class \code{"\linkS4class{SnpCallSetPlus}"}, directly.
46
-Class \code{"\linkS4class{SnpQSet}"}, by class "SnpCallSetPlus", distance 2.
47
-Class \code{"\linkS4class{SnpCallSet}"}, by class "SnpCallSetPlus", distance 2.
48
-Class \code{"\linkS4class{QuantificationSet}"}, by class "SnpCallSetPlus", distance 3.
49
-Class \code{"\linkS4class{SnpLevelSet}"}, by class "SnpCallSetPlus", distance 3.
50
-Class \code{"\linkS4class{eSet}"}, by class "SnpCallSetPlus", distance 4.
51
-Class \code{"\linkS4class{VersionedBiobase}"}, by class "SnpCallSetPlus", distance 5.
52
-Class \code{"\linkS4class{Versioned}"}, by class "SnpCallSetPlus", distance 6.
53
-}
54
-\section{Methods}{
55
-  \describe{
56
-    \item{CA}{\code{signature(object = "CNSet")}: ... }
57
-    \item{CA<-}{\code{signature(object = "CNSet", value = "matrix")}: ... }
58
-    \item{CB}{\code{signature(object = "CNSet")}: ... }
59
-    \item{CB<-}{\code{signature(object = "CNSet", value = "matrix")}: ... }
60
-    \item{computeHmm}{\code{signature(object = "CNSet")}: ... }
61
-    \item{copyNumber}{\code{signature(object = "CNSet")}: ... }
62
-    \item{ellipse}{\code{signature(x = "CNSet")}: ... }
63
-	 }
64
-}
65
-\author{R. Scharpf}
66
-\examples{
67
-showClass("CNSet")
68
-}
69
-\keyword{classes}
70 0
deleted file mode 100644
... ...
@@ -1,41 +0,0 @@
1
-\name{SnpSet-methods}
2
-\docType{methods}
3
-%\alias{calls}
4
-\alias{calls,SnpSet-method}
5
-\alias{calls<-,SnpSet,matrix-method}
6
-%\alias{confs}
7
-\alias{confs,SnpSet-method}
8
-%\alias{confs<-}
9
-\alias{confs<-,SnpSet,matrix-method}
10
-\title{Accessors for Calls and Confidences on a SnpSet object}
11
-
12
-\description{
13
-  \code{calls} returns the genotype calls. CRLMM stores genotype calls
14
-  as integers (1 - AA; 2 - AB; 3 - BB).
15
-
16
-  \code{confs} returns the confidences associated to the genotype
17
-  calls. The current implementation of CRLMM stores the confidences as
18
-  integers by using the transformation:
19
-
20
-  conf = round(-1000*log2(1-p)),
21
-
22
-  where 'p' is the posterior probability of the call.
23
-}
24
-
25
-\section{Methods}{
26
-\describe{
27
-  \item{\code{initialize(SnpSet)}:}{Object instantiation, used by
28
-    \code{new}; not be be called directly by the user.}  
29
-  \item{\code{calls(object)}}{accessor for genotype calls}
30
-  \item{\code{confs(object)}}{accessor for crlmm genotype confidence scores}  
31
-}}
32
-
33
-\examples{
34
-  theCalls <- matrix(sample(1:3, 20, rep=TRUE), nc=2)
35
-  p <- matrix(runif(20), nc=2)
36
-  theConfs <- round(-1000*log2(1-p))
37
-  obj <- new("SnpSet", call=theCalls, callProbability=theConfs)
38
-  calls(obj)
39
-  confs(obj)
40
-}
41
-\keyword{manip}
42 0
new file mode 100644
... ...
@@ -0,0 +1,298 @@
1
+\name{cnOptions}
2
+\alias{cnOptions}
3
+\title{
4
+	Options for copy number estimation
5
+}
6
+\description{
7
+	This function returns all the user-modifiable arguments to the crlmm copy number function.
8
+}
9
+\usage{
10
+    cnOptions(outdir = "./", cdfName, crlmmFile = "snpsetObject.rda", intensityFile = "normalizedIntensities.rda", rgFile = "rgFile.rda", save.it = TRUE, save.cnset = TRUE, load.it = TRUE, splitByChr = TRUE, MIN.OBS = 3, MIN.SAMPLES = 10, batch = NULL, DF.PRIOR = 50, bias.adj = FALSE, prior.prob = rep(1/4, 4), SNRmin = 4, chromosome = 1:24, seed = 123, verbose = TRUE, GT.CONF.THR = 0.99, PHI.THR = 2^6, nHOM.THR = 5, MIN.NU = 2^3, MIN.PHI = 2^3, THR.NU.PHI = TRUE, thresholdCopynumber = TRUE, unlink = TRUE, ...)
11
+}
12
+\arguments{
13
+  \item{outdir}{
14
+	Path to store output from genotyping / copy number algorithms
15
+}
16
+  \item{cdfName}{
17
+        Character string indicating array-type.  See crlmm:::validCdfNames().
18
+}
19
+  \item{crlmmFile}{
20
+
21
+        When \code{save.it} is \code{TRUE}, the output from the crlmm
22
+        genotyping will be saved to \code{<outdir>/<crlmmFile>}.  This
23
+        object should not be directly loaded by the user.  When
24
+        \code{load.it} is \code{TRUE}, the function
25
+        \code{crlmmCopynumber} will load this object from
26
+        \code{outdir} without rerunning the quantile normalization and
27
+        genotyping steps. 
28
+
29
+}
30
+
31
+  \item{intensityFile}{
32
+
33
+        When \code{save.it} is \code{TRUE}, the allele summaries of
34
+        the quantile normalized intensities at polymorphic loci will
35
+        be saved to \code{<outdir>/<intensityFile>}.  This object is
36
+        not intended to be loaded directly by the user.  When
37
+        \code{load.it} is \code{TRUE}, the function
38
+        \code{crlmmCopynumber} will load this object from
39
+        \code{outdir} without rerunning the quantile normalization and
40
+        genotyping steps.
41
+       
42
+}
43
+
44
+  \item{rgFile}{
45
+	
46
+	For Affymetrix platforms, \code{rgFile} is ignored.  When
47
+        \code{save.it} is \code{TRUE}, the R and G summaries of the
48
+        quantile normalized intensities for illumina platforms is
49
+        saved to \code{<outdir>/<rgFile>}.  This object is not
50
+        intended to be loaded directly by the user.  When
51
+        \code{load.it} is \code{TRUE}, the function
52
+        \code{crlmmCopynumber} will load this object from
53
+        \code{outdir} without rerunning the quantile normalization and
54
+        genotyping steps.
55
+
56
+}
57
+
58
+  \item{save.it}{
59
+  
60
+  When \code{TRUE}, intermediate files containing quantile-normalized
61
+  intensities and genotyping results are saved to \code{outdir}.
62
+  Saving these objects can save time if the copy number estimation is
63
+  repeated.  Intermediate files for the Affymetrix platform are
64
+  \code{intensityFile} and \code{crlmmFile}; for the Illumina
65
+  platform, \code{rgfile} is an intermediate file.
66
+  
67
+}
68
+
69
+  \item{save.cnset}{
70
+
71
+  This argument is ignored for the processing of more than one
72
+  chromosome (e.g., \code{length(chromosome) > 1}).  If \code{TRUE},
73
+  results from the copy number estimation for each chromosome is saved
74
+  to a file.  The format of the saved file is:
75
+
76
+  \code{<outdir>/cnSet_<chromosome>.rda}
77
+
78
+}
79
+
80
+  \item{load.it}{
81
+
82
+  If \code{TRUE}, intermediate files are loaded from \code{outdir}.
83
+  See \code{rgFile}, \code{intensityFile}, \code{crlmmFile}.
84
+
85
+}
86
+
87
+  \item{splitByChr}{
88
+  
89
+  If \code{TRUE}, results are saved for each chromosome.  In general,
90
+  this should always be \code{TRUE}.
91
+
92
+}
93
+
94
+  \item{MIN.OBS}{
95
+
96
+  For genotypes with fewer than \code{MIN.OBS}, the within-genotype
97
+  median is imputed from the observed within-genotype mediants at that
98
+  loci.  The parameters used in the regression are estimated from
99
+  polymorphic loci where the genotypes AA, AB, and BB each have a
100
+  frequency greater than \code{MIN.OBS}.
101
+  
102
+}
103
+
104
+  \item{MIN.SAMPLES}{
105
+
106
+  The minimum number of samples required for each batch. Batches
107
+  wither fewer than \code{MIN.SAMPLES} are skipped.
108
+
109
+}
110
+
111
+  \item{batch}{
112
+
113
+  Character string or factor denoting the batch for each file to be
114
+  processed.  The length of this argument should be the same as the
115
+  number of files. The batch covariate is a surrogate for experimental
116
+  conditions that change over calendar time. Typically, batch can be
117
+  denoted by the 96 well chemistry plate or the month / year.  
118
+
119
+}
120
+
121
+  \item{DF.PRIOR}{
122
+
123
+  The 2 x 2 covariance matrix of the background and signal variances
124
+  is estimated from the data at each locus.  This matrix is then
125
+  smoothed towards a common matrix estimated from all of the loci.
126
+  DF.PRIOR controls the amount of smoothing towards the common matrix,
127
+  with higher values corresponding to greater smoothing.  Currently,
128
+  DF.PRIOR is not estimated from the data.  Future versions may
129
+  estimate DF.PRIOR empirically.
130
+
131
+}
132
+
133
+  \item{bias.adj}{
134
+
135
+  If \code{TRUE}, initial estimates of the linear model are updated
136
+  after excluding samples that have a low posterior probability of
137
+  normal copy number.  Excluding samples that have a low posterior
138
+  probability can be helpful at loci in which a substantial fraction
139
+  of the samples have a copy number alteration.  For additional
140
+  information, see Scharpf et al., 2009.
141
+
142
+}
143
+
144
+  \item{prior.prob}{
145
+
146
+  A numerical vector providing prior probabilities for copy number
147
+  states corresponding to homozygous deletion, hemizygous deletion,
148
+  normal copy number, and amplification, respectively.
149
+
150
+}
151
+
152
+  \item{SNRmin}{
153
+
154
+  The signal to noise ratio (SNR) estimated during the CRLMM
155
+  genotyping is a summary measure of sample quality based on the
156
+  separation of the genotype clusters.  Smaller values of the SNR
157
+  correspond to samples of lower quality.  Samples are excluded from
158
+  the copy number estimation step if SNR values are less than this
159
+  value. A SNR less than 5 for the Affymetrix platform generally
160
+  corresponds to low quality.  For the Illumina platform, we have
161
+  observed samples with poor quality with SNR < 32.  More specific
162
+  recommendations for Illumina are being evaluated.
163
+
164
+}
165
+
166
+  \item{chromosome}{
167
+
168
+  The chromosome(s) to estimate copy number.  Valid entries are 1-23,
169
+  where 23 corresponds to chromosome X.  Quantile normalization and
170
+  genotyping are performed for all SNPs and nonpolymorphic features,
171
+  irrespective of the \code{chromosome} argument.
172
+
173
+}
174
+
175
+  \item{seed}{
176
+
177
+  Seed for random number generation (integer).  Used only for reproducibility.
178
+
179
+}
180
+
181
+  \item{verbose}{
182
+
183
+  If \code{TRUE}, verbose output.
184
+
185
+}
186
+
187
+  \item{GT.CONF.THR}{
188
+
189
+    Confidence threshold for genotype calls (0, 1).  Calls with
190
+    confidence scores below this theshold are not used to estimate the
191
+    within-genotype medians.
192
+
193
+}
194
+
195
+  \item{PHI.THR}{
196
+
197
+    SNPs with slopes (phi values) below this value are flagged.
198
+    Flagged SNPs are not used in a regression to impute background and
199
+    slope coefficients at nonpolymorphic loci.
200
+
201
+}
202
+
203
+  \item{nHOM.THR}{
204
+
205
+  If fewer than \code{nHOM.THR} homozygous genotypes (AA or BB) are
206
+    observed, the SNPs is flagged.  Flagged SNPs are not used in a
207
+    regression to impute background and slope coefficients at
208
+    nonpolymorphic loci.
209
+
210
+}
211
+
212
+  \item{MIN.NU}{
213
+
214
+  Minimum value for the estimate of background in the linear model.
215
+  Negative values are permissible from the estimation, but not
216
+  plausible.  Ignored if \code{THR.NU.PHI} is \code{FALSE}.
217
+
218
+}
219
+
220
+  \item{MIN.PHI}{
221
+
222
+  Minimum value for the estimate of signal (phi) in the linear model.
223
+  Negative values are permissible from the estimation, but not
224
+  plausible.   Ignored if \code{THR.NU.PHI} is \code{FALSE}.
225
+
226
+}
227
+
228
+  \item{THR.NU.PHI}{
229
+
230
+  If \code{THR.NU.PHI} is \code{FALSE}, \code{MIN.NU} and
231
+  \code{MIN.PHI} are ignored.
232
+
233
+}
234
+
235
+  \item{thresholdCopynumber}{
236
+
237
+  If \code{TRUE}, allele-specific number estimates are truncated.
238
+  Values less than 0.05 are assigned the value 0.05; values exceeding
239
+  5 are assigned the value 5.  
240
+
241
+}
242
+
243
+  \item{unlink}{
244
+
245
+  Whether to remove intermediate files storing the quantile normalized
246
+  intensities.  Ignored if \code{save.it} is \code{FALSE}.
247
+
248
+}
249
+
250
+  \item{\dots}{
251
+
252
+  Additional arguments are passed to \code{readIdatFiles} (Illumina
253
+  platforms only).
254
+
255
+}
256
+}
257
+\details{
258
+
259
+	The minimum required arguments when calling this function are
260
+	\code{cdfName} and \code{batch}.  The user will generally want
261
+	to assign a valid path to \code{outdir} that specifies where
262
+	the intermediate files and processed data are saved.
263
+	
264
+}
265
+
266
+\value{
267
+
268
+	In general, nothing is returned and all results are saved to
269
+	\code{outdir}.
270
+
271
+}
272
+
273
+
274
+\references{
275
+
276
+	RB Scharpf, I Ruczinski, B Carvalho, B Doan, A Chakravarti,
277
+	and R Irizarry (2009), A multilevel model to address batch
278
+	effects in copy number estimation using SNP arrays (Technical
279
+	Report).
280
+
281
+}
282
+
283
+\author{  R. Scharpf}
284
+
285
+
286
+\seealso{
287
+}
288
+\examples{
289
+require(hapmapsnp6)
290
+path <- system.file("celFiles", package="hapmapsnp6")
291
+celfiles <- list.celfiles(path)
292
+## the different populations were run in different batches.  For the
293
+##  files in this package, the batch is indicated by the 13th character
294
+##  in the string
295
+batch <- substr(celfiles, 13, 13)
296
+}
297
+\keyword{manip}
298
+
0 299
deleted file mode 100644
... ...
@@ -1,20 +0,0 @@
1
-\name{methods-eSet}
2
-%%\alias{chromosome}
3
-\alias{chromosome,eSet-method}
4
-%%\alias{position}
5
-\alias{position,eSet-method}
6
-\title{Methods for eSet derivatives}
7
-
8
-\description{
9
-  \code{chromosome} return the chromosome number for SNP and NP probes.
10
-  \code{position} return the physical position of the SNP or the
11
-  physical position for the first index of the NP probe.
12
-}
13
-
14
-\section{Methods}{
15
-  \describe{
16
-    \item{\code{chromosome(object)}:}{ accessor for chromosome}
17
-    \item{\code{position(object)}:}{accessor for genomic position}
18
-  }
19
-}
20
-\keyword{manip}
21 0
\ No newline at end of file
22 1
deleted file mode 100644
... ...
@@ -1,45 +0,0 @@
1
-\name{list.celfiles}
2
-\alias{list.celfiles}
3
-
4
-\title{List CEL files.}
5
-\description{
6
-  Function used to get a list of CEL files.
7
-}
8
-\usage{
9
-list.celfiles(...)
10
-}
11
-
12
-\arguments{
13
-  \item{\dots}{Same arguments of \code{\link{list.files}}}
14
-}
15
-\details{
16
-  For the moment, this function returns only uncompressed CEL files (ie,
17
-  no CEL.gz)
18
-}
19
-\value{
20
-  Character vector with filenames.
21
-}
22
-
23
-\note{
24
-  Quite often users want to use this function to pass filenames to other
25
-  methods. In this situations, it is safer to use the argument 'full.names=TRUE'.
26
-}
27
-\seealso{\code{\link{list.files}}}
28
-\examples{
29
-if (require(hapmapsnp5)){
30
-  path <- system.file("celFiles", package="hapmapsnp5")
31
-
32
-  ## only the filenames
33
-  list.celfiles(path)
34
-
35
-  ## the filenames with full path...
36
-  ## very useful when genotyping samples not in the working directory
37
-  list.celfiles(path, full.names=TRUE)
38
-}else{
39
-  ## this won't return anything
40
-  ## if in the working directory there isn't any CEL
41
-  list.celfiles(getwd())
42
-}
43
-}
44
-\keyword{IO}
45
-\keyword{utilities}
... ...
@@ -50,7 +50,7 @@ snprma(filenames, mixtureSampleSize = 10^5, fitMixture = FALSE, eps = 0.1, verbo
50 50
   \item{cdfName}{Name of the CDF}
51 51
 }
52 52
 \examples{
53
-if (require(genomewidesnp5Crlmm) & require(hapmapsnp5)){
53
+if (require(genomewidesnp5Crlmm) & require(hapmapsnp5) & require(oligoClasses)){
54 54
   path <- system.file("celFiles", package="hapmapsnp5")
55 55
 
56 56
   ## the filenames with full path...