@@ -290,3 +290,18 @@ is decoded and scanned

 * updated addFeatureAnnotation so that CHR arg not required

 * fixed bug in nonpolymorphic function -- function checks whether any regions are pseudoautosomal 

 * fixed bug in list2locusset function (no longer assigns genomewidesnp6 as the default annotation)

+

+2009-10-16 R. Scharpf - committed version 1.3.23

+

+ * 100*CA, 100*CB are stored in the CA, CB assayDataElements.  The

+   replacement method automatically converts a matrix on the copy

+   number scale to an integer representation by scaling by a factor of

+   100. The accessor methods for CA, CB automatically divide by 100 to

+   return estimates back on the original copy number scale.

+

+ * more debugging of readIdatFiles and illumina copy number vignette.

+

+ * updated illumina_copynumber vignette

+

+

+

@@ -1,13 +1,13 @@

 Package: crlmm

 Type: Package

 Title: Genotype Calling (CRLMM) and Copy Number Analysis tool for Affymetrix SNP 5.0 and 6.0 and Illumina arrays.

-Version: 1.3.22

+Version: 1.3.23

 Date: 2009-09-29

 Author: Rafael A Irizarry, Benilton S Carvalho <bcarvalh@jhsph.edu>, Robert Scharpf <rscharpf@jhsph.edu>, Matt Ritchie <mritchie@wehi.EDU.AU>

 Maintainer: Benilton S Carvalho <bcarvalh@jhsph.edu>, Robert Scharpf <rscharpf@jhsph.edu>, Matt Ritchie <mritchie@wehi.EDU.AU>

 Description: Faster implementation of CRLMM specific to SNP 5.0 and 6.0 arrays, as well as a copy number tool specific to 6.0.

 License: Artistic-2.0

-Depends: methods, Biobase (>= 2.5.5)

+Depends: methods, Biobase (>= 2.5.5), R (>= 2.10.0)

 Imports: affyio, preprocessCore, utils, stats, genefilter, splines, mvtnorm, oligoClasses, ellipse, methods, SNPchip

 Suggests: hapmapsnp5, hapmapsnp6, genomewidesnp5Crlmm (>= 1.0.2),genomewidesnp6Crlmm (>= 1.0.2), snpMatrix, VanillaICE (>= 1.7.8), GGdata

 Collate: AllClasses.R

@@ -217,15 +217,16 @@ harmonizeDimnamesTo <- function(object1, object2){

 }

 crlmmWrapper <- function(filenames,

-			 cdfName="genomewidesnp6",

+			 cdfName,

 			 load.it=FALSE,

 			 save.it=FALSE,

 			 splitByChr=TRUE,

 			 crlmmFile,

 			 intensityFile,

 			 rgFile,

-			 platform=c("affymetrix", "illumina")[1],

 			 ...){

+	if(missing(cdfName)) stop("cdfName is missing -- a valid cdfName is required.  See crlmm:::validCdfNames()")  

+	platform <- whichPlatform(cdfName)

 	if(!(platform %in% c("affymetrix", "illumina"))){

 		stop("Only 'affymetrix' and 'illumina' platforms are supported at this time.")

 	} else {

@@ -234,7 +235,11 @@ crlmmWrapper <- function(filenames,

 	if(missing(intensityFile)) stop("must specify 'intensityFile'.")

 	if(missing(crlmmFile)) stop("must specify 'crlmmFile'.")

 	if(platform == "illumina"){

-		if(missing(rgFile)) stop("must specify 'rgFile'.")

+		if(missing(rgFile)){

+			##stop("must specify 'rgFile'.")

+			rgFile <- file.path(dirname(crlmmFile), "rgFile.rda")

+			message("rgFile not specified.  Using ", rgFile)

+		}

 		if(!load.it){

 			RG <- readIdatFiles(...)

 			if(save.it) save(RG, file=rgFile)

@@ -264,7 +269,7 @@ crlmmWrapper <- function(filenames,

 	}

 	if(load.it){

 		if(!file.exists(crlmmFile)){

-			message("load.it is TRUE, but 'crlmmFile' does not exist.  Rerunning the genotype calling algorithm") 

+			message("load.it is TRUE, but ", crlmmFile, " does not exist.  Rerunning the genotype calling algorithm") 

 			load.it <- FALSE

 		}

 	}

@@ -603,7 +608,11 @@ computeCopynumber <- function(object,

 	if(ncol(object) < 10)

 		stop("Must have at least 10 samples in each batch to estimate model parameters....preferably closer to 90 samples per batch")

 	##require(oligoClasses)

-	if(missing(CHR)) stop("Must specify CHR")

+	if(missing(CHR)){

+		if(length(unique(chromosome(object)))==1){

+			CHR <- unique(chromosome(object))

+		} else  stop("Must specify CHR")

+	}

 	if(CHR == 24) stop("Nothing available yet for chromosome Y")

 	if(missing(batch)) {

 		message("'batch' missing.  Assuming all samples in the CrlmmSetList object were processed together in the same batch.")

@@ -803,7 +812,6 @@ list2locusSet <- function(envir, ABset, NPset, CHR, cdfName){

 	dimnames(CA) <- list(envir[["snps"]], envir[["sns"]])	

 	CB <- envir[["CB"]]

 	dimnames(CB) <- dimnames(CA)

-

 	##NP copy number

 	if(!is.null(NPset)){

 		CT <- envir[["CT"]]

@@ -907,15 +915,19 @@ list2locusSet <- function(envir, ABset, NPset, CHR, cdfName){

 	rownames(fd) <- rownames(CA)[I]

 	fD <- new("AnnotatedDataFrame",

 		  data=fd,

-		  varMetadata=data.frame(labelDescription=colnames(fd)))	

+		  varMetadata=data.frame(labelDescription=colnames(fd)))

+	placeholder <- matrix(NA, nrow(CA), ncol(CA))

+	dimnames(placeholder) <- dimnames(CA)

 	assayData <- assayDataNew("lockedEnvironment",

-				  CA=CA[I, ],

-				  CB=CB[I, ])

+				  CA=placeholder[I, ],

+				  CB=placeholder[I, ])

 	cnset <- new("CopyNumberSet",

 		      assayData=assayData,

 		      featureData=fD,

 		      phenoData=phenodata,

 		      annotation=cdfName)

+	CA(cnset) <- CA  ##replacement method converts to integer

+	CB(cnset) <- CB  ##replacement method converts to integer

 	cnset <- thresholdCopyNumberSet(cnset)

 	return(cnset)

 }

@@ -929,10 +941,11 @@ thresholdCopyNumberSet <- function(object){

 	ca[ca > 5] <- 5

 	cb[cb < 0.05] <- 0.05

 	cb[cb > 5] <- 5

-	ca <- matrix(as.integer(ca*100), nrow(ca), ncol(ca))

-	cb <- matrix(as.integer(cb*100), nrow(cb), ncol(cb))

-	rownames(ca) <- rownames(cb) <- featureNames(object)

-	colnames(ca) <- colnames(cb) <- sampleNames(object)

+##	ca <- matrix(as.integer(ca*100), nrow(ca), ncol(ca))

+##	cb <- matrix(as.integer(cb*100), nrow(cb), ncol(cb))

+##	rownames(ca) <- rownames(cb) <- featureNames(object)

+##	colnames(ca) <- colnames(cb) <- sampleNames(object)

+	##replacement method multiplies by 100 and converts to integer

 	CA(object) <- ca

 	CB(object) <- cb

 	return(object)

@@ -1143,8 +1156,8 @@ nonpolymorphic <- function(plateIndex, NP, envir, CONF.THR=0.99, DF.PRIOR=50, pk

 		}

 		CT1 <- 1/phi1*(NP[, gender=="male"]-nu1)

 		CT2 <- 1/phi2*(NP[, gender=="female"]-nu2)

-		CT1 <- matrix(as.integer(100*CT1), nrow(CT1), ncol(CT1))

-		CT2 <- matrix(as.integer(100*CT2), nrow(CT2), ncol(CT2))

+		##CT1 <- matrix(as.integer(100*CT1), nrow(CT1), ncol(CT1))

+		##CT2 <- matrix(as.integer(100*CT2), nrow(CT2), ncol(CT2))

 		CT <- envir[["CT"]]

 		CT[, plate==uplate[p] & gender=="male"] <- CT1

 		CT[, plate==uplate[p] & gender=="female"] <- CT2

@@ -1163,7 +1176,7 @@ nonpolymorphic <- function(plateIndex, NP, envir, CONF.THR=0.99, DF.PRIOR=50, pk

 		normalNP <- envir[["normalNP"]]

 		normalNP <- normalNP[, plate==uplate[p]]

 		mus <- rowMedians(NP * normalNP, na.rm=TRUE)

-		crosshyb <- median(muA) - median(mus)

+		crosshyb <- max(median(muA) - median(mus), 0)

 		X <- cbind(1, log2(mus+crosshyb))

 		##X <- cbind(1, log2(mus))

 		logPhiT <- X %*% betahat

@@ -1171,7 +1184,8 @@ nonpolymorphic <- function(plateIndex, NP, envir, CONF.THR=0.99, DF.PRIOR=50, pk

 		nuT[, p] <- mus - 2*phiT[, p]

 		T <- 1/phiT[, p]*(NP - nuT[, p])

 		CT <- envir[["CT"]]

-		CT[, plate==uplate[p]] <- matrix(as.integer(100*T), nrow(T), ncol(T))

+##		CT[, plate==uplate[p]] <- matrix(as.integer(100*T), nrow(T), ncol(T))

+		CT[, plate==uplate[p]] <- matrix(T, nrow(T), ncol(T))

 		##Variance for prediction region

 		##sig2T[, plate==uplate[[p]]] <- rowMAD(log2(NP*normalNP), na.rm=TRUE)^2

@@ -1615,11 +1629,11 @@ polymorphic <- function(plateIndex, A, B, envir){

 		tmp <- (B-nuB - phistar*A + phistar*nuA)/phiB

 		copyB <- tmp/(1-phistar*phiAx/phiB)

 		copyA <- (A-nuA-phiAx*copyB)/phiA

-		CB[, plate==uplate[p]] <- matrix(as.integer(100*copyB), nrow(copyB), ncol(copyB))

-		CA[, plate==uplate[p]] <- matrix(as.integer(100*copyA), nrow(copyA), ncol(copyA))

+		CB[, plate==uplate[p]] <- copyB

+		CA[, plate==uplate[p]] <- copyA

 	} else{

-		CA[, plate==uplate[p]] <- matrix(as.integer(100*1/phiA*(A-nuA)), nrow(A), ncol(A))

-		CB[, plate==uplate[p]] <- matrix(as.integer(100*1/phiB*(B-nuB)), nrow(A), ncol(A))

+		CA[, plate==uplate[p]] <- matrix((1/phiA*(A-nuA)), nrow(A), ncol(A))

+		CB[, plate==uplate[p]] <- matrix((1/phiB*(B-nuB)), nrow(B), ncol(B))

 	}

 	assign("CA", CA, envir)

 	assign("CB", CB, envir)

@@ -1979,7 +1993,8 @@ setMethod("update", "character", function(object, ...){

 		if(!"chromosome" %in% fvarLabels(crlmmSetList[[1]])){

 			featureData(crlmmSetList[[1]]) <- addFeatureAnnotation(crlmmSetList)

 		} 

-		CHR <- unique(chromosome(crlmmSetList[[1]]))		

+		CHR <- unique(chromosome(crlmmSetList[[1]]))

+		if(length(CHR) > 1) stop("More than one chromosome in the object. This method requires one chromosome at a time.")		

 		if(CHR==24){

 			message("skipping chromosome 24")

 			next()

@@ -21,7 +21,7 @@ readIdatFiles <- function(sampleSheet=NULL,

 	       pd = new("AnnotatedDataFrame", data = data.frame(Sample_ID=arrayNames))

        }	

        if(!is.null(sampleSheet)) { # get array info from Illumina's sample sheet

-	       if(!is.null(arrayNames)){

+	       if(is.null(arrayNames)){

 		       ##arrayNames=NULL

 		       if(!is.null(arrayInfoColNames$barcode) && (arrayInfoColNames$barcode %in% colnames(sampleSheet))) {

 			       barcode = sampleSheet[,arrayInfoColNames$barcode]

@@ -41,7 +41,7 @@ readIdatFiles <- function(sampleSheet=NULL,

 		       }

 	       }

 	       pd = new("AnnotatedDataFrame", data = sampleSheet)

-	       sampleNames(pd) <- arrayNames

+	       sampleNames(pd) <- basename(arrayNames)

        }

        narrays = length(arrayNames)

        grnfiles = paste(arrayNames, fileExt$green, sep=sep)

@@ -6,10 +6,19 @@ setValidity("CopyNumberSet", function(object) {

 setMethod("CA", "CopyNumberSet", function(object, ...) assayData(object)[["CA"]]/100)

 setMethod("CB", "CopyNumberSet", function(object, ...) assayData(object)[["CB"]]/100)

-setReplaceMethod("CA", signature(object="CopyNumberSet", value="matrix"),

-		 function(object, value) assayDataElementReplace(object, "CA", value*100))

-setReplaceMethod("CB", signature(object="CopyNumberSet", value="matrix"),

-		 function(object, value) assayDataElementReplace(object, "CB", value*100))

+setReplaceMethod("CA", signature(object="CopyNumberSet", value="matrix"), function(object, value){

+	dns <- dimnames(value)

+	value <- matrix(as.integer(value*100), nrow(value), ncol(value))

+	dimnames(value) <- dns

+	assayDataElementReplace(object, "CA", value)

+})

+

+setReplaceMethod("CB", signature(object="CopyNumberSet", value="matrix"), function(object, value){

+	dns <- dimnames(value)	

+	value <- matrix(as.integer(value*100), nrow(value), ncol(value))

+	dimnames(value) <- dns	

+	assayDataElementReplace(object, "CB", value)

+})

@@ -73,6 +73,12 @@ setMethod("addFeatureAnnotation", "CrlmmSetList", function(object, ...){

 	position <- c(position.snp, position.np)

 	chrom <- c(chr.snp, chr.np)

+

+	##We may not have annotation for all of the snps

+	if(!all(featureNames(object) %in% names(position))){

+		message("Dropping loci for which physical position  is not available.")

+		object <- object[featureNames(object) %in% names(position), ]

+	}

 	ix <- match(featureNames(object), names(position))

 	position <- position[ix]

 	chrom <- chrom[ix]

@@ -200,8 +206,19 @@ setMethod("splitByChromosome", "CrlmmSetList", function(object, cdfName, outdir)

 		save(crlmmSetList, file=file.path(outdir, paste("crlmmSetList_", CHR, ".rda", sep="")))

 	}

 })

+

 setMethod("update", "CrlmmSetList", function(object, ...){

-	computeCopynumber(object, ...)

+	if(length(crlmmSetList) == 3){

+		message("copy number object already present. Nothing to do.")

+		return()

+	}

+	CHR <- unique(chromosome(crlmmSetList[[1]]))

+	if(length(CHR) > 1) stop("More than one chromosome in the object. This method requires one chromosome at a time.")

+	if(CHR == 24){

+		message("A solution for chromosome 24 is not yet available.")

+		return()

+	}	

+	computeCopynumber(object, CHR=CHR, ...)

 })

 setReplaceMethod("CA", signature(object="CrlmmSetList", value="matrix"),

@@ -47,6 +47,11 @@ library(crlmm)

 crlmm:::validCdfNames()

 @ 

+<<Rsge, eval=TRUE, echo=FALSE>>=

+##If multiple cluster nodes are available

+library(Rsge)

+@ 

+

 \paragraph{Preprocess and genotype.}

 Specify the coordinates of Affymetrix cel files and where to put

@@ -56,11 +61,19 @@ number estimation.

 <<celfiles>>=

 myPath <- "/thumper/ctsa/snpmicroarray/hapmap/raw/affy/1m"

 celFiles <- list.celfiles(myPath, full.names=TRUE, pattern=".CEL")

-#outdir <- "/thumper/ctsa/snpmicroarray/rs/data/hapmap/1m/affy"

-outdir <- "/thumper/ctsa/snpmicroarray/rs/data/hapmap/1m/tmp"

-dir.create(outdir)

+outdir <- "/thumper/ctsa/snpmicroarray/rs/data/hapmap/1m/affy"

+#outdir <- "/thumper/ctsa/snpmicroarray/rs/data/hapmap/1m/tmp"

+if(!file.exists(outdir)) dir.create(outdir)

 @ 

+Split the list of files by plate.  For this HapMap data, the different

+ancestries were run on separate plates.

+

+<<plate>>=

+plate <- substr(basename(celFiles), 13, 13)

+table(plate)

+celFiles <- split(celFiles, plate)

+@ 

 \noindent The next code chunk quantile normalizes the samples to a

@@ -68,21 +81,76 @@ target reference distribution and uses the crlmm algorithm to genotype.

 See \cite{Carvalho2007a} for details regarding the crlmm genotyping

 algorithm.

-<<preprocessAndGenotype, eval=FALSE>>=

-crlmmWrapper(celFiles,##[1:15], 

-	     cdfName="genomewidesnp6",

-	     save.it=TRUE, 

-	     load.it=FALSE,

-	     intensityFile=file.path(outdir, "normalizedIntensities.rda"),

-	     crlmmFile=file.path(outdir, "snpsetObject.rda"),

-	     platform="affymetrix")

-@ 

-

-As a result of the above wrapper, the following R objects are now created:

+<<preprocessAndGenotype, eval=TRUE, echo=TRUE>>=

+for(i in seq(along=celFiles)){

+	platedir <- file.path(outdir, names(celFiles)[i])

+	if(!file.exists(platedir)) dir.create(platedir)

+	crlmmWrapper(celFiles[[i]],##[1:15], 

+		     cdfName="genomewidesnp6",

+		     save.it=FALSE, 

+		     load.it=TRUE,

+		     intensityFile=file.path(platedir, "normalizedIntensities.rda"),

+		     crlmmFile=file.path(platedir, "snpsetObject.rda"))

+}

+q("no")

+@ 

+

+<<preprocessAndGenotype, eval=FALSE, echo=FALSE>>=

+for(i in seq(along=celFiles)){

+	platedir <- file.path(outdir, names(celFiles)[i])

+	if(!file.exists(platedir)) dir.create(platedir)

+	sge.submit(crlmmWrapper,

+		   filenames=celFiles[[i]],##[1:15], 

+		   cdfName="genomewidesnp6",

+		   save.it=FALSE, 

+		   load.it=TRUE,

+		   intensityFile=file.path(platedir, "normalizedIntensities.rda"),

+		   crlmmFile=file.path(platedir, "snpsetObject.rda"),

+		   global.savelist=c("platedir"),

+		   packages="crlmm")

+}

+##platedirs <- sapply(names(celFiles), function(x) file.path(outdir, x))

+##intensityFiles <- file.path(platedirs, "normalizedIntensities.rda")

+####crlmmFiles <- file.path(platedirs, "snpsetObject.rda")

+####myF <- function(celFiles, intensityFile, crlmmFile, ...){

+####	crlmmWrapper(celFiles, intensityFile=intensityFile,

+####		     crlmmFile=crlmmFile,

+##sge.apply(celFiles, 

+##	  crlmmWrapper,

+##	  cdfName="genomewidesnp6",

+##	  save.it=FALSE, 

+##	  load.it=TRUE,

+##	  intensityFile=,

+##	  crlmmFile=file.path(platedir, "snpsetObject.rda"))

+q("no")

+@ 

+

+As a result of the above wrapper, the following R objects are now

+created for each plate:

+

+<<crlmmSetListObjects>>=

+fns <- list.files(platedir, pattern="crlmmSetList", full.name=TRUE)

+@ 

+

+To estimate copy number, simply run the update function on this vector

+of filenames.  If multiple nodes are available for computation, one

+could easily split this job.

+

+<<copynumber, echo=FALSE, eval=TRUE>>=

+for(i in seq(along=celFiles)){

+	fns <- list.files(file.path(outdir, names(celFiles[i])), 

+			  pattern="crlmmSetList", full.names=TRUE)

+	update(fns)

+}

+@

-<<>>=

-list.files(outdir)[grep("crlmmSetList", list.files(outdir))]

-@ 

+<<copynumber.sge, echo=FALSE, eval=TRUE>>=

+for(i in seq(along=celFiles)){

+	fns <- list.files(file.path(outdir, names(celFiles[i])), 

+			  pattern="crlmmSetList", full.names=TRUE)

+	sge.parLapply(X=fns, FUN=update, packages="crlmm", cluster=FALSE)

+}

+@

 \paragraph{Locus- and allele-specific estimates of copy number.}

@@ -92,9 +160,7 @@ Load the object for chromosome 22 and compute copy number:

 CHR <- 22

 if(!exists("crlmmSetList")) load(file.path(outdir, paste("crlmmSetList_", CHR, ".rda", sep="")))

 show(crlmmSetList)

-if(length(crlmmSetList)==2){

-	crlmmSetList <- update(crlmmSetList, CHR=CHR)

-}

+crlmmSetList <- update(crlmmSetList)

 show(crlmmSetList)

 @ 

@@ -501,10 +567,10 @@ for(i in indices[[j]]){

     intensities.  Each panel displays a single SNP.}

 \end{figure}

-TODO: Plot the prediction regions for total copy number 2 and 3 for the

-first plate. Plotting symbols are the genotype calls (1=AA, 2=AB, 3=BB);

-light grey points are from other plates. One could also add the

-prediction regions for 0-4 copies, but it gets crowded.

+% TODO: Plot the prediction regions for total copy number 2 and 3 for

+% the first plate. Plotting symbols are the genotype calls (1=AA, 2=AB,

+% 3=BB); light grey points are from other plates. One could also add the

+% prediction regions for 0-4 copies, but it gets crowded.

 <<predictionRegion, fig=TRUE, width=8, height=8, include=FALSE, eval=FALSE, echo=FALSE>>=

@@ -572,8 +638,10 @@ the allele-specific copy number.  SNP-specific parameters are estimated

 only from samples with a suitable confidence score for the diallelic

 genotype calls.  The default confidence threshold is 0.99, but can be

 adjusted by passing the argument \Robject{CONF.THR} to the

-\Robject{update} method.  TODO: illustrate this approach with boxplots

-of the A and B intensities stratified by genotype.

+\Robject{update} method.  

+

+% TODO: illustrate this approach with boxplots of the A and B

+% intensities stratified by genotype.

 <<linearModel, fig=TRUE,include=FALSE, eval=FALSE, echo=FALSE>>=

 cols <- brewer.pal(7, "Accent")[c(1, 4, 7)]

@@ -28,50 +28,70 @@ options(prompt="R> ")

 <<>>=

 library(crlmm)

+crlmm:::validCdfNames()

+outdir <- "/thumper/ctsa/snpmicroarray/rs/data/hapmap/illumina/HumanCNV370-Duo"

+datadir <- "/thumper/ctsa/snpmicroarray/illumina/IDATS/370k"

+cdfName <- "human370v1c"

 @ 

-<<echo=FALSE, eval=TRUE>>=

-setwd("/thumper/ctsa/snpmicroarray/illumina/IDATS/370k")

-@ 

-

-<<readIdat>>=

-samplesheet5 = read.csv("HumanHap370Duo_Sample_Map.csv", header=TRUE, as.is=TRUE)[-c(28:46,61:75,78:79),]

-if(!exists("RG")){

-	RG <- readIdatFiles(samplesheet5, arrayInfoColNames=list(barcode=NULL, position="SentrixPosition"), saveDate=TRUE)

-}

+<<samplesToProcess>>=

+samplesheet = read.csv(file.path(datadir, "HumanHap370Duo_Sample_Map.csv"), header=TRUE, as.is=TRUE)

+## Excluding a few of the problem samples

+samplesheet <- samplesheet[-c(28:46,61:75,78:79), ]

+arrayNames <- file.path(datadir, unique(samplesheet[, "SentrixPosition"]))

+##Check that files exist

+grnfiles = all(file.exists(paste(arrayNames, ".Grn.dat", sep="")))

+redfiles = all(file.exists(paste(arrayNames, ".Red.dat", sep="")))

 @ 

 Alternatively, arguments to the readIdatFiles can be passed through the

 {\tt ...} argument of the \R{} function \Robject{crlmmWrapper}.

 <<wrapper>>=

-crlmmWrapper(sampleSheet=samplesheet5,

+crlmmWrapper(sampleSheet=samplesheet,

+	     arrayNames=arrayNames,

 	     arrayInfoColNames=list(barcode=NULL, position="SentrixPosition"), 

 	     saveDate=TRUE,

 	     cdfName=cdfName,

 	     load.it=FALSE,

-	     save.it=TRUE,

-	     intensityFile=file.path(platedir, "normalizedIntensities.rda"),

-	     crlmmFile=file.path(platedir, "snpsetObject.rda"),

-	     rgFile=file.path(platedir, "rgFile.rda"),

-	     platform="illumina")

+	     save.it=FALSE,

+	     intensityFile=file.path(outdir, "normalizedIntensities.rda"),

+	     crlmmFile=file.path(outdir, "snpsetObject.rda"),

+	     rgFile=file.path(outdir, "rgFile.rda"))

 @ 

-Samples with low signal to noise ratios tend to have a lot of variation

-in the point estimates of copy number.  One may want to exclude these

-samples.

+This creates a \Robject{crlmmSetList} object in the \Robject{outdir}

+directory.  The first element of this object contains the

+quantile-normalized A and B intensities.  The second element in the list

+contains the crlmm genotype calls.  

 <<SNR>>=

-hist(crlmmOut$SNR)

+CHR <- 1

+filename <- paste(outdir, "/crlmmSetList_", CHR, ".rda", sep="")

+load(filename)

+hist(crlmmSetList[[1]]$SNR)

 @ 

-Run update on the CrlmmSetList object to obtain copy number. Estimate

+Run update on the CrlmmSetList object to obtain copy number estimates. Estimate

 copy number for chromosome 1.

-<<>>=

-CHR <- 1

-filename <- paste(platedir, "/CrlmmSetList_", CHR, ".rda", sep="")

+<<estimateCn>>=

+update(filename)

+##Or equivalently,

+##crlmmSetList <- computeCopynumber(crlmmSetList)

+##save(crlmmSetList, file=file.path(outdir, paste("crlmmSetList_", CHR, ".rda", sep="")))

+@ 

+

+Samples with low signal to noise ratios tend to have a lot of variation

+in the point estimates of copy number.  One may want to exclude these

+samples, or smooth after filtering outliers.  Here we load the

+crlmmSetList object.  See the copynumber.Rnw vignette for example plots.

+

+<<crlmmSetList>>=

+##reload the updated object

+load(filename)

+cn <- copyNumber(crlmmSetList)

 @ 

 \end{document}

new file mode 100644

Binary files /dev/null and b/inst/scripts/illumina_copynumber.pdf differ

@@ -115,7 +115,7 @@

     \item{"CA"}{\code{signature(object = "CrlmmSetList")}: extracts the

       copy number for allele A at polymorphic loci. For nonpolymorphic

-      probes, CA returns the total copy number.}

+      probes, CA returns the total copy number. }

     \item{"CB"}{\code{signature(object = "CrlmmSetList")}: extracts the

       copy number for allele B at polymorphic loci. For nonpolymorphic

@@ -9,9 +9,8 @@

   confidence scores are assigned using the crlmm algorithm.

 }

 \usage{

-crlmmWrapper(filenames, cdfName = "genomewidesnp6", load.it=FALSE,

-save.it=FALSE, splitByChr=TRUE, crlmmFile, intensityFile, rgFile,

-platform=c("affymetrix", "illumina")[1], ...)

+crlmmWrapper(filenames, cdfName, load.it=FALSE,

+save.it=FALSE, splitByChr=TRUE, crlmmFile, intensityFile, rgFile, ...)

 }

 \arguments{

   \item{filenames}{

@@ -19,6 +18,7 @@ platform=c("affymetrix", "illumina")[1], ...)

   }

   \item{cdfName}{

     Annotation package (currently, only 'genomewidesnp6' is supported).

+    See \code{crlmm:::validCdfNames()}.

   }

   \item{load.it}{ Logical.  If TRUE, the wrapper function will try to

     load the 'intensityFile' and 'crlmmFile'. When FALSE, the quantile

@@ -41,8 +41,6 @@ platform=c("affymetrix", "illumina")[1], ...)

   \item{rgFile}{Character string.  Name of file (including the

     complete path) containing the RG intensities (only applicable for

     illumina platform).  }

-  \item{platform}{Character string.  Currently only \code{affymetrix}

-    and \code{illumina} are supported.}

   \item{\dots}{

     additional arguments to function \code{readIdatFiles} (illumina

     platform only)