@@ -264,3 +264,11 @@ is decoded and scanned

 * overhaul of vignette

 * added update() method for copy number

 * added documentation and error checks for crlmmWrapper

+

+2009-08-14 R. Scharpf - committed version 1.3.19

+

+* changes to crlmmWrapper (trying to make this handle the illumina platform as well...still needs testing)

+  - removed crlmmIlluminaWrapper.Rd, and commented crlmmIlluminaWrapper function

+* labeled figures / displayed output of code chunks in the copy number vignette

+* added bibliography for copy number vignette.  Added file inst/doc/refs.bib

+* added boxplot method 

@@ -1,7 +1,7 @@

 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.18

+Version: 1.3.19

 Date: 2009-07-22

 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>

@@ -28,7 +28,7 @@ importFrom(Biobase, assayDataElement, assayDataElementNames,

            assayDataElementReplace, assayDataNew, classVersion, validMsg)

 importFrom(graphics, abline, axis, layout, legend, mtext, par, plot,

-           polygon, rect, segments, text, points)

+           polygon, rect, segments, text, points, boxplot)

 importFrom(stats, update)

@@ -76,7 +76,7 @@ export(batch,

        copyNumber,

        crlmm,

        crlmmIllumina,

-       crlmmIlluminaWrapper,

+##       crlmmIlluminaWrapper,

        crlmmWrapper,

 ##       ellipse,

 ##       computeEmission,

@@ -215,69 +215,76 @@ harmonizeDimnamesTo <- function(object1, object2){

 	return(object1)

 }

-crlmmIlluminaWrapper <- function(outdir="./",

-				 cdfName,

-				 save.intermediate=FALSE,

-				 splitByChr=TRUE,

-				 intensityFile,

-				 ...){  ##additional arguments to readIdatFiles

-	stopifnot(basename(intensityFile) == "res.rda")

-	if(!file.exists(outdir)) stop(outdir, " does not exist.")

-	if(!isValidCdfName(cdfName, platform="illumina")) stop(cdfName, " not supported.")

-	if(file.exists(file.path(outdir, "RG.rda"))) {

-		message("Loading RG.rda...")

-		load(file.path(outdir, "RG.rda"))

-	} else {

-		message("Reading Idat files...")		

-		RG <- readIdatFiles(...)

-		##J <- match(c("1_A", "3_A", "5_A", "7_A"), sampleNames(RG))

-		##RG <- RG[, -J]

-		if(save.intermediate) save(RG, file=file.path(outdir, "RG.rda"))  ##935M for 91 samples...better not to save this

-	}	

-	if(!file.exists(intensityFile)){

-		message("Quantile normalization / genotyping...")				

-		crlmmOut <- crlmmIllumina(RG=RG,

-					  cdfName=cdfName,

-                                          sns=sampleNames(RG),

-                                          returnParams=TRUE,

-                                          save.it=TRUE,

-                                          intensityFile=intensityFile)

-		if(save.intermediate) save(crlmmOut, file=file.path(outdir, "crlmmOut.rda"))

-	} else{

-		load(file.path(outdir, "crlmmOut.rda"))

-	}

-	message("Loading ", intensityFile, "...")		

-	load(intensityFile)

-	if(exists("cnAB")){

-		np.A <- as.integer(cnAB$A)

-		np.B <- as.integer(cnAB$B)

-		np <- ifelse(np.A > np.B, np.A, np.B)

-		np <- matrix(np, nrow(cnAB$A), ncol(cnAB$A))

-		rownames(np) <- cnAB$gns

-		colnames(np) <- cnAB$sns

-		cnAB$NP <- np

-	}

-	sampleNames(crlmmOut) <- res$sns	

-	if(exists("cnAB")){

-		ABset <- combineIntensities(get("res"), cnAB, cdfName=cdfName)

-	} else{

-		ABset <- combineIntensities(get("res"), NULL, cdfName=cdfName)

-	}

-	##protocolData(ABset)[["ScanDate"]] <- as.character(pData(RG)$ScanDate)

-	crlmmResult <- harmonizeSnpSet(crlmmOut, ABset, cdfName)

-	stopifnot(all.equal(dimnames(crlmmOut), dimnames(ABset)))

-	crlmmList <- list(ABset,

-			  crlmmResult)

-	crlmmList <- as(crlmmList, "CrlmmSetList")

-	if(splitByChr){

-		message("Saving by chromosome")

-		splitByChromosome(crlmmList, cdfName=cdfName, outdir=outdir)

-	} else{

-		message("Saving crlmmList object to ", outdir)

-		save(crlmmList, file=file.path(outdir, "crlmmList.rda"))

-	}

-	message("CrlmmSetList objects saved to ", outdir)

-}

+##crlmmIlluminaWrapper <- function(filenames,

+##				 cdfName,

+##				 load.it=FALSE,

+##				 save.it=FALSE,

+##				 splitByChr=TRUE,

+##				 crlmmFile,

+##				 intensityFile, ...){

+####				 outdir="./",

+####				 cdfName,

+####				 save.intermediate=FALSE,

+####				 splitByChr=TRUE,

+####				 intensityFile,

+####				 ...){  ##additional arguments to readIdatFiles

+####	stopifnot(basename(intensityFile) == "res.rda")

+##	if(!file.exists(outdir)) stop(outdir, " does not exist.")

+##	if(!isValidCdfName(cdfName, platform="illumina")) stop(cdfName, " not supported.")

+##	if(file.exists(file.path(outdir, "RG.rda"))) {

+##		message("Loading RG.rda...")

+##		load(file.path(outdir, "RG.rda"))

+##	} else {

+##		message("Reading Idat files...")		

+##		RG <- readIdatFiles(...)

+##		##J <- match(c("1_A", "3_A", "5_A", "7_A"), sampleNames(RG))

+##		##RG <- RG[, -J]

+##		if(save.intermediate) save(RG, file=file.path(outdir, "RG.rda"))  ##935M for 91 samples...better not to save this

+##	}	

+##	if(!file.exists(intensityFile)){

+##		message("Quantile normalization / genotyping...")				

+##		crlmmOut <- crlmmIllumina(RG=RG,

+##					  cdfName=cdfName,

+##                                          sns=sampleNames(RG),

+##                                          returnParams=TRUE,

+##                                          save.it=TRUE,

+##                                          intensityFile=intensityFile)

+##		if(save.intermediate) save(crlmmOut, file=file.path(outdir, "crlmmOut.rda"))

+##	} else{

+##		load(file.path(outdir, "crlmmOut.rda"))

+##	}

+##	message("Loading ", intensityFile, "...")		

+##	load(intensityFile)

+##	if(exists("cnAB")){

+##		np.A <- as.integer(cnAB$A)

+##		np.B <- as.integer(cnAB$B)

+##		np <- ifelse(np.A > np.B, np.A, np.B)

+##		np <- matrix(np, nrow(cnAB$A), ncol(cnAB$A))

+##		rownames(np) <- cnAB$gns

+##		colnames(np) <- cnAB$sns

+##		cnAB$NP <- np

+##	}

+##	sampleNames(crlmmOut) <- res$sns	

+##	if(exists("cnAB")){

+##		ABset <- combineIntensities(get("res"), cnAB, cdfName=cdfName)

+##	} else{

+##		ABset <- combineIntensities(get("res"), NULL, cdfName=cdfName)

+##	}

+##	##protocolData(ABset)[["ScanDate"]] <- as.character(pData(RG)$ScanDate)

+##	crlmmResult <- harmonizeSnpSet(crlmmOut, ABset, cdfName)

+##	stopifnot(all.equal(dimnames(crlmmOut), dimnames(ABset)))

+##	crlmmList <- list(ABset,

+##			  crlmmResult)

+##	crlmmList <- as(crlmmList, "CrlmmSetList")

+##	if(splitByChr){

+##		message("Saving by chromosome")

+##		splitByChromosome(crlmmList, cdfName=cdfName, outdir=outdir)

+##	} else{

+##		message("Saving crlmmList object to ", outdir)

+##		save(crlmmList, file=file.path(outdir, "crlmmList.rda"))

+##	}

+##	message("CrlmmSetList objects saved to ", outdir)

+##}

 ##quantileNormalize1m <- function(cel.files,

 ##				outdir,

@@ -393,7 +400,29 @@ crlmmWrapper <- function(filenames,

 			 save.it=FALSE,

 			 splitByChr=TRUE,

 			 crlmmFile,

-			 intensityFile, ...){

+			 intensityFile,

+			 rgFile,

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

+			 ...){

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

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

+	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(load.it & !file.exists(rgFile)){

+			message("load.it is TRUE, bug rgFile not present.  Attempting to read the idatFiles.")

+			RG <- readIdatFiles(...)

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

+		}

+		if(load.it & file.exists(rgFile)){

+			message("Loading RG file")

+			load(rgFile)

+			RG <- get("RG")

+		}

+	}

+	if(!(file.exists(dirname(crlmmFile)))) stop(dirname(crlmmFile), " does not exist.")

+	if(!(file.exists(dirname(intensityFile)))) stop(dirname(intensityFile), " does not exist.")

 	outfiles <- file.path(dirname(crlmmFile), paste("crlmmSetList_", 1:24, ".rda", sep=""))

 	if(load.it & all(file.exists(outfiles))){

 		load(outfiles[1])

@@ -404,38 +433,55 @@ crlmmWrapper <- function(filenames,

 			return("load.it is TRUE and 'crlmmSetList_<CHR>.rda' objects found. Nothing to do...")

 		}

 	}

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

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

 	if(load.it){

 		if(!file.exists(crlmmFile)){

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

 			load.it <- FALSE

 		}

 	}

-	if(!file.exists(outdir)) stop(outdir, " does not exist.")

-	if(!isValidCdfName(cdfName, platform="affymetrix"))

-		stop(cdfName, " is not a valid entry.  See crlmm:::validCdfNames(platform='affymetrix')")

-	if(!file.exists(crlmmFile) | !load.it){

-		crlmmResult <- crlmm(filenames=filenames,

-				     cdfName=cdfName,

-				     save.it=TRUE,

-				     load.it=load.it,

-				     intensityFile=intensityFile, ...)

-		message("Quantile normalizing the copy number probes...")		

-		cnrmaResult <- cnrma(filenames=filenames, cdfName=cdfName)

-		if(save.it){

-			message("Saving crlmmResult and cnrmaResult to", crlmmFile)

-			save(crlmmResult, cnrmaResult, file=crlmmFile)

+	if(!isValidCdfName(cdfName, platform=platform))

+		stop(cdfName, " is not a valid entry.  See crlmm:::validCdfNames(platform)")

+	if(platform == "affymetrix"){

+		if(!file.exists(crlmmFile) | !load.it){

+			crlmmResult <- crlmm(filenames=filenames,

+					     cdfName=cdfName,

+					     save.it=TRUE,

+					     load.it=load.it,

+					     intensityFile=intensityFile)

+			message("Quantile normalizing the copy number probes...")		

+			cnrmaResult <- cnrma(filenames=filenames, cdfName=cdfName)

+			if(save.it){

+				message("Saving crlmmResult and cnrmaResult to", crlmmFile)

+				save(crlmmResult, cnrmaResult, file=crlmmFile)

+			}

+		} else {

+			message("Loading ", crlmmFile, "...")

+			load(crlmmFile)

+			crlmmResult <- get("crlmmResult")

+			cnrmaResult <- get("cnrmaResult")

+		}

+	}

+	if(platform == "illumina"){

+		if(!file.exists(crlmmFile) | !load.it){		

+			crlmmOut <- crlmmIllumina(RG=RG,

+						  cdfName=cdfName,

+						  sns=sampleNames(RG),

+						  returnParams=TRUE,

+						  save.it=TRUE,

+						  intensityFile=intensityFile)

+			if(save.it) save(crlmmOut, file=crlmmFile)

+		} else {

+			message("Loading ", crlmmFile, "...")

+			load(crlmmFile)

+			crlmmResult <- get("crlmmResult")

+			if(exists("cnAB")){

+				cnrmaResult <- get("cnAB")

+			} else cnrmaResult <- NULL

 		}

-	} else {

-		message("Loading ", crlmmFile, "...")

-		load(crlmmFile)

-		crlmmResult <- get("crlmmResult")

-		cnrmaResult <- get("cnrmaResult")

 	}

 	load(intensityFile)

 	ABset <- combineIntensities(get("res"), cnrmaResult, cdfName)

-	protocolData(ABset)[["ScanDate"]] <- as.character(celDates(filenames))	

+	if(platform=="affymetrix") protocolData(ABset)[["ScanDate"]] <- as.character(celDates(filenames))	

 	crlmmResult <- harmonizeSnpSet(crlmmResult, ABset, cdfName)

 	stopifnot(all.equal(dimnames(crlmmResult), dimnames(ABset)))

 	crlmmSetList <- as(list(ABset, crlmmResult), "CrlmmSetList")

@@ -611,9 +657,7 @@ instantiateObjects <- function(calls, conf, NP, plate, envir,

 	envir[["normalNP"]] <- normalNP

 }

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

-	computeCopynumber(object, ...)

-})

+

 computeCopynumber <- function(object,

 			      CHR,

@@ -621,6 +665,7 @@ computeCopynumber <- function(object,

 			      batch,

 			      SNRmin=5,

 			      cdfName, ...){

+	if(class(object) != "CrlmmSetList") stop("object must be of class ClrmmSetList")

 	if(missing(cdfName))

 		cdfName <- annotation(object)

 	if(!isValidCdfName(cdfName, platform="affymetrix")) stop(cdfName, " not supported.")	

@@ -23,43 +23,48 @@ setMethod("[", "CrlmmSetList", function(x, i, j, ..., drop = FALSE){

 	as(x, "CrlmmSetList")	

 })

-##setReplaceMethod("[[", "CrlmmSetList", function(x, i, j, ..., value) {

-##	browser()

-##	##otherwise infinite recursion

-##	x <- as(x, "list")

-##	x[[i]] <- value

-##	x <- as(x, "CrlmmSetList")

-##	x <- .harmonizeDimnames(x)

-##	stopifnot(identical(featureNames(x[[i]]), featureNames(x[[1]])))

-##	return(x)

-##})

+setMethod("$", "CrlmmSetList", function(x, name) {

+	##if(!(name %in% .parameterNames()[output(x) != 0])){

+	if(length(x) != 3){

+		stop("'$' operature reserved for accessing parameter names in CopyNumberSet object.  CrlmmSetList must be of length 3")

+	}

+	j <- grep(name, fvarLabels(x[[3]]))	

+	if(length(j) < 1)

+		stop(name, " not in fvarLabels of CopyNumberSet object")

+	if(length(j) > 1){

+		warning("Multiple instances of ", name, " in fvarLabels.  Using the first instance")

+		j <- j[1]

+	}

+	param <- fData(x[[3]])[, j]

+	param

+})

+setMethod(".harmonizeDimnames", "CrlmmSetList", function(object){

+	i <- length(object)

+	while(i > 1){

+		object[[i-1]] <- harmonizeDimnamesTo(object[[i-1]], object[[i]])

+		i <- i-1

+	}

+	object

+})

 setMethod("A", "CrlmmSetList", function(object) A(object[[1]]))

-

 setMethod("annotation", "CrlmmSetList", function(object) annotation(object[[1]]))

-

 setMethod("B", "CrlmmSetList", function(object) B(object[[1]]))

-

+setMethod("batch", "CrlmmSetList", function(object) batch(object[[3]]))

 setMethod("CA", "CrlmmSetList", function(object, ...) CA(object[[3]], ...))

 setMethod("CB", "CrlmmSetList", function(object, ...) CB(object[[3]], ...))

-

 setMethod("calls", "CrlmmSetList", function(object) calls(object[[2]]))

+setMethod("chromosome", "CrlmmSetList", function(object) chromosome(object[[3]]))

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

 	match(cnNames(object[[1]], annotation(object)), featureNames(object))

 })

-

-setMethod("copyNumber", "CrlmmSetList", function(object) copyNumber(object[[3]]))

-

 setMethod("combine", signature=signature(x="CrlmmSetList", y="CrlmmSetList"),

           function(x, y, ...){

 		  x.abset <- x[[1]]

 		  y.abset <- y[[1]]

-

 		  x.snpset <- x[[2]]

 		  y.snpset <- y[[2]]

-

 		  abset <- combine(x.abset, y.abset)

-

 		  ##we have hijacked the featureData slot to store parameters.  Biobase will not allow combining our 'feature' data.

 		  warning("the featureData is not easily combined...  removing the featureData")

 		  ##fd1 <- featureData(x.snpset)

@@ -71,13 +76,10 @@ setMethod("combine", signature=signature(x="CrlmmSetList", y="CrlmmSetList"),

 		  merged <- as(merged, "CrlmmSetList")

 		  merged

 	  })

-		  

-

-

-##setMethod("fData", "CrlmmSetList", function(object) featureNames(object[[1]]))

-

+setMethod("confs", "CrlmmSetList", function(object) confs(object[[2]]))

+setMethod("copyNumber", "CrlmmSetList", function(object) copyNumber(object[[3]]))

+setMethod("dims", "CrlmmSetList", function(object) sapply(object, dim))

 setMethod("featureNames", "CrlmmSetList", function(object) featureNames(object[[1]]))

-

 setMethod("ncol", signature(x="CrlmmSetList"), function(x) ncol(x[[1]]))

 setMethod("nrow", signature(x="CrlmmSetList"), function(x) nrow(x[[1]]))

 setMethod("plot", signature(x="CrlmmSetList"),

@@ -86,18 +88,15 @@ setMethod("plot", signature(x="CrlmmSetList"),

 		  B <- log2(B(x))

 		  plot(A, B, ...)

 	  })

-

 setMethod("points", signature(x="CrlmmSetList"),

 	  function(x, y, ...){

 		  A <- log2(A(x))

 		  B <- log2(B(x))

 		  points(A, B, ...)

 	  })

-

+setMethod("position", "CrlmmSetList", function(object) position(object[[3]]))

 setMethod("sampleNames", "CrlmmSetList", function(object) sampleNames(object[[1]]))

-

 setMethod("show", "CrlmmSetList", function(object){

-	##for(i in seq(along=object)) show(object[[i]])

 	cat("\n Elements in CrlmmSetList object: \n")

 	cat("\n")

 	for(i in 1:length(object)){

@@ -106,43 +105,7 @@ setMethod("show", "CrlmmSetList", function(object){

 		cat("Dimensions: ", dim(object[[i]]))		

 		cat("\n \n")

 	}

-##	cat("\n")

-##	cat("Dimensions:\n")

-##	print(dims(object))

-})

-

-setMethod("chromosome", "CrlmmSetList", function(object) chromosome(object[[3]]))

-setMethod("position", "CrlmmSetList", function(object) position(object[[3]]))

-setMethod("confs", "CrlmmSetList", function(object) confs(object[[2]]))

-

-

-setMethod(".harmonizeDimnames", "CrlmmSetList", function(object){

-	i <- length(object)

-	while(i > 1){

-		object[[i-1]] <- harmonizeDimnamesTo(object[[i-1]], object[[i]])

-		i <- i-1

-	}

-	object

-})

-setMethod("dims", "CrlmmSetList", function(object) sapply(object, dim))

-setMethod("batch", "CrlmmSetList", function(object) batch(object[[3]]))

-setMethod("$", "CrlmmSetList", function(x, name) {

-	##if(!(name %in% .parameterNames()[output(x) != 0])){

-	if(length(x) != 3){

-		stop("'$' operature reserved for accessing parameter names in CopyNumberSet object.  CrlmmSetList must be of length 3")

-	}

-	j <- grep(name, fvarLabels(x[[3]]))	

-	if(length(j) < 1)

-		stop(name, " not in fvarLabels of CopyNumberSet object")

-	if(length(j) > 1){

-		warning("Multiple instances of ", name, " in fvarLabels.  Using the first instance")

-		j <- j[1]

-	}

-	param <- fData(x[[3]])[, j]

-	param

 })

-

-	

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

 	match(snpNames(object[[1]], annotation(object)), featureNames(object))

 })

@@ -163,5 +126,78 @@ 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, ...)

+})

+setMethod("boxplot", "CrlmmSetList", function(x, ...){

+##boxplot.CrlmmSetList <- function(x, ...){

+	if(length(x) != 3) stop("elements of list should be of class ABset, SnpSet, and CopyNumberSet, respectively.")

+	genotypes <- calls(x)-1

+	A1 <- A(x)

+	B1 <- B(x)

+	Alist <- split(A1, genotypes)

+	Alist <- rev(Alist)

+	Blist <- split(B1, genotypes)

+	ylim <- range(unlist(Alist))

+	boxplot(Alist, xaxt="n", ylab=expression(I[A]), 

+		cex.axis=0.6,

+		xlab="",

+		ylim=range(unlist(Alist), na.rm=TRUE),	

+		border="grey50", xaxs="i", at=0:2, xlim=c(-0.5, 2.5),

+		cex.main=0.9, xaxt="n",

+		yaxt="n",

+		col=cols)

+	axis(2, at=pretty(ylim), cex.axis=0.8)

+	axis(1, at=0:2, labels=rev(c("2A (AA genotype)", "1A (AB genotype)", "0A (BB genotype)")), cex.axis=0.7)

+	##extracts nuA for first batch

+	suppressWarnings(nuA <- x$nuA)

+	segments(-1, nuA, 

+		 0, nuA, lty=2, col="blue")

+	suppressWarnings(phiA <- x$phiA)

+	##phiA <- fData(x)[["phiA_A"]]

+	segments(0, nuA,

+		 2.5, nuA+2.5*phiA, lwd=2, col="blue")

+	axis(2, at=nuA, labels=expression(hat(nu[A])), cex.axis=0.9)

+	text(0, ylim[1], labels=paste("n =", length(Alist[[1]])), cex=0.8)

+	text(1, ylim[1], labels=paste("n =", length(Alist[[2]])), cex=0.8)

+	text(2, ylim[1], labels=paste("n =", length(Alist[[3]])), cex=0.8)

+##	segments(0.5, nuA+0.5*phiA,

+##		 1, pretty(ylim, n=5)[2], lty=2, col="blue")

+##	segments(1, pretty(ylim, n=5)[2],

+##		 1.2, pretty(ylim)[2], lty=2, col="blue")

+##	text(1.25, pretty(ylim)[2], pos=4, 

+##	     labels=expression(hat(nu[A]) + c[A]*hat(phi[A])))

+	legend("topleft", fill=rev(cols), legend=c("AA", "AB", "BB"))##, title="diallelic genotypes")	

+

+	ylim <- range(unlist(Blist))

+	boxplot(Blist, xaxt="n", ylab=expression(I[B]), 

+		##xlab="diallelic genotypes", 

+		cex.axis=0.6, 

+		ylim=range(unlist(Blist), na.rm=TRUE),

+		border="grey50", xaxs="i", 

+		at=0:2, xlim=c(-0.5, 2.5),

+		cex.main=0.9, yaxt="n", 

+		col=rev(cols))

+	axis(2, at=pretty(ylim), cex.axis=0.8)

+	axis(1, at=0:2, labels=c("0B (AA genotype)", "1B (AB genotype)", "2B (BB genotype)"), cex.axis=0.7)

+	##nuB <- fData(x)[["nuB_A"]]

+	suppressWarnings(nuB <- x$nuB)

+	segments(-1, nuB, 

+		 0, nuB, lty=2, col="blue")

+	##phiB <- fData(x[i,])[["phiB_A"]]

+	suppressWarnings(phiB <- x$phiB)

+	segments(0, nuB,

+		 2.5, nuB+2.5*phiB, lwd=2, col="blue")

+	axis(2, at=nuB, labels=expression(hat(nu[B])), cex.axis=0.9)

+	text(0, ylim[1], labels=paste("n =", length(Blist[[1]])), cex=0.8)

+	text(1, ylim[1], labels=paste("n =", length(Blist[[2]])), cex=0.8)

+	text(2, ylim[1], labels=paste("n =", length(Blist[[3]])), cex=0.8)

+##	segments(0.5, nuB+0.5*phiB,

+##		 1, pretty(ylim,n=5)[2], lty=2, col="blue")

+##	segments(1, pretty(ylim, n=5)[2],

+##		 1.2, pretty(ylim,n=5)[2], lty=2, col="blue")

+##	text(1.25, pretty(ylim,n=5)[2], pos=4, labels=expression(hat(nu[B]) + c[B]*hat(phi[B])))	

+	mtext(featureNames(x)[i], 3, outer=TRUE, line=1)	

+})

@@ -4,6 +4,7 @@

 %\VignettePackage{crlmm}

 \documentclass{article}

 \usepackage{graphicx}

+\usepackage[authoryear,round,numbers]{natbib}

 \newcommand{\Rfunction}[1]{{\texttt{#1}}}

 \newcommand{\Rmethod}[1]{{\texttt{#1}}}

 \newcommand{\Rcode}[1]{{\texttt{#1}}}

@@ -32,7 +33,9 @@ options(prompt="R> ")

 \begin{abstract}

   This vignette estimates copy number for HapMap samples on the

-  Affymetrix 6.0 platform.

+  Affymetrix 6.0 platform.  See \citep{Scharpf2009} for the working

+  paper.

+  

 \end{abstract}

 \section{Simple usage}

@@ -142,7 +145,32 @@ str(brks)

 @ 

 \paragraph{Circular binary segmentation}

-TO DO.

+

+<<setup>>=

+library(DNAcopy)

+library(genefilter)

+ratioset <- crlmmSetList[[3]]

+meds <- rowMedians(copyNumber(ratioset), na.rm=TRUE)

+ratios <- log2(copyNumber(ratioset)) - log2(meds)

+ratioset <- new("SnpCopyNumberSet",

+		copyNumber=ratios,

+		phenoData=phenoData(ratioset),

+		featureData=featureData(ratioset),

+		annotation=annotation(ratioset))

+@ 

+

+The following code chunk (not evaluated) takes a while to run:

+

+<<cbs, eval=FALSE>>=

+CNA.object <- CNA(copyNumber(ratioset), 

+		  chromosome(ratioset),

+		  position(ratioset),

+		  data.type="logratio",

+		  sampleid=sampleNames(ratioset))

+smoothed.CNA.object <- smooth.CNA(CNA.object)

+segment.cna <- segment(smoothed.CNA.object, verbose=1)

+save(segment.cna, file=file.path(outdir, paste("segment.cna_", CHR, ".rda", sep="")))

+@ 

 \section{Accessors}

@@ -301,7 +329,12 @@ argument to \Robject{TRUE}:

 update(crlmmSetList, CHR=22, bias.adj=TRUE)

 @ 

-Calculate the frequency of amplifications and deletions at each locus.

+The following code chunk calculates the frequency of amplifications and

+deletions at each locus. Shaded regions above the zero line in Figure

+\ref{fig:hmm_hapmap} display the frequency of amplifications, whereas

+shaded regions below the zero line graphically display the frequency of

+hemizygous or homozygous deletions.

+

 <<hmm_hapmap, fig=TRUE, include=FALSE, width=8, height=7>>=

 require(SNPchip)

 library(RColorBrewer)

@@ -351,14 +384,17 @@ gc()

 \begin{figure}

   \centering

   \includegraphics[width=\textwidth]{copynumber-hmm_hapmap}

-  \caption{}

+  \caption{\label{fig:hmm_hapmap} The frequency of amplifications in the

+    hapmap samples is displayed above the zero line.  The frequency of

+    hemizygous or homozygous deletions are displayed below the zero

+    line.}

 \end{figure}

+\clearpage

 \paragraph{One sample at a time: locus-level estimates}

-Plot physical position versus copy number for the first sample.  Recall

-that the copy number estimates were multiplied by 100 and stored as an

-integer.

+Figure \ref{fig:oneSample} plots physical position (horizontal axis)

+versus copy number (vertical axis) for the first sample.  

 <<oneSample, fig=TRUE, width=8, height=4, include=FALSE>>=

 par(las=1)

@@ -377,17 +413,55 @@ plotCytoband(22, new=FALSE, cytoband.ycoords=c(5.8, 6.0), label.cytoband=FALSE)

 \begin{figure}

   \includegraphics[width=0.9\textwidth]{copynumber-oneSample}

-  \caption{Total copy number (y-axis) for chromosome 22 plotted

-    against physical position (x-axis) for one sample.  Estimates at

-    nonpolymorphic loci are plotted in light blue.}

+  \caption{\label{fig:oneSample} Total copy number (y-axis) for

+    chromosome 22 plotted against physical position (x-axis) for one

+    sample.  Estimates at nonpolymorphic loci are plotted in light

+    blue.}

+\end{figure}

+

+\noindent The following code chunk plots the locus-level copy number

+estimates and overlays the predictions from the hidden markov model.

+

+<<overlayHmmPredictions, fig=TRUE, include=FALSE>>=

+ask <- FALSE

+op <- par(mfrow=c(3, 1), las=1, mar=c(1, 4, 1, 1), oma=c(3, 1, 1, 1), ask=ask)

+##Put fit on the copy number scale

+cns <- copyNumber(crlmmSetList)

+cnState <- hmmPredictions - as.integer(1)

+xlim <- c(10*1e6, max(position(crlmmSetList)))

+cols <- brewer.pal(8, "Dark2")[1:4]

+for(j in 1:3){

+	plot(position(crlmmSetList), cnState[, j], pch=".", col=cols[2], xaxt="n", 

+	     ylab="copy number", xlab="Physical position (Mb)", type="s", lwd=2,

+	     ylim=c(0,6), xlim=xlim)

+	points(position(crlmmSetList), cns[, j], pch=".", col=cols[3])

+	lines(position(crlmmSetList), cnState[,j], lwd=2, col=cols[2])

+	axis(1, at=pretty(position(crlmmSetList)), 

+	     labels=pretty(position(crlmmSetList))/1e6)

+	abline(h=c(1,3), lty=2, col=cols[1])	

+	legend("topright", bty="n", legend=sampleNames(crlmmSetList)[j])

+	legend("topleft", lty=1, col=cols[2], legend="copy number state",

+	       bty="n", lwd=2)

+	plotCytoband(CHR, cytoband.ycoords=c(5, 5.2), new=FALSE,

+		     label.cytoband=FALSE, xlim=xlim)

+}

+par(op)

+@ 

+

+\begin{figure}

+  \includegraphics[width=\textwidth]{copynumber-overlayHmmPredictions}

+  \caption{\label{fig:overlayHmmPredictions} Total copy number (y-axis)

+    for chromosome 22 plotted against physical position (x-axis) for

+    three samples.  Estimates at nonpolymorphic loci are plotted in

+    light blue. }

 \end{figure}

+\clearpage

 \paragraph{One SNP at a time}

-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.

+Scatterplots of the A and B allele intensities (log-scale) can be useful

+for assessing the diallelic genotype calls.  The following code chunk is

+displayed in Figure \ref{fig:genotypeCalls}.

 <<genotypeCalls, fig=TRUE, width=8, height=8, include=FALSE>>=

 xlim <- ylim <- c(6.5,13)

@@ -397,7 +471,7 @@ cex <- 0.6

 par(mfrow=c(3,3), las=1, pty="s", ask=FALSE, mar=c(2, 2, 2, 2), oma=c(2, 2, 1, 1))

 ##plot 9 at a time

 indices <- split(snpIndex(crlmmSetList), rep(1:length(snpIndex(crlmmSetList)), each=9, length.out=length(snpIndex(crlmmSetList))))

-par(ask=FALSE)

+##par(ask=FALSE)

 j <- 1

 for(i in indices[[j]]){

 	gt <- calls(crlmmSetList)[i, ]

@@ -409,9 +483,22 @@ for(i in indices[[j]]){

 	mtext("A", 1, outer=TRUE, line=1)

 	mtext("B", 2, outer=TRUE, line=1)	

 }

-@ 

+@

+

+\begin{figure}

+  \centering

+  \includegraphics[width=0.8\textwidth]{copynumber-genotypeCalls}

+  \caption{\label{fig:genotypeCalls}  Scatterplots of A versus B

+    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.

+

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

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

 require(RColorBrewer)

 library(ellipse)

 greens <- brewer.pal(9, "Greens")

@@ -470,14 +557,45 @@ j <- 1

 %\end{figure}

 \section{Details for the copy number estimation procedure}

-TO DO.

+

+We assume a linear relationship between the normalized intensities and

+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.

+

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

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

+i <- match("SNP_A-8608180", featureNames(crlmmSetList))

+B <- batch(crlmmSetList) == unique(batch(crlmmSetList))[1]

+gt.conf <- confs(crlmmSetList[i, B])

+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))

+crlmm:::boxplot(crlmmSetList[i, B][[1]])

+par(op)

+@ 

+

+%\begin{figure}

+%  \centering

+%  \includegraphics[width=0.8\textwidth]{copynumber-linearModel}

+%

+%  \caption{\label{fig:linearModel} The SNP-specific intercept and slope

+%    for the first batch of samples.  These plots can sometimes be

+%    misleading as the genotypes estimated with low confidence do not

+%    influence the regression line.}

+%\end{figure}

 \section{Session information}

 <<sessionInfo, results=tex>>=

 toLatex(sessionInfo())

 @ 

+\section*{References}

-

+%\begin{bibliography}

+  \bibliographystyle{plain}

+  \bibliography{refs}

+%\end{bibliography}

 \end{document}

Binary files a/inst/scripts/copynumber.pdf and b/inst/scripts/copynumber.pdf differ

new file mode 100644

@@ -0,0 +1,8 @@

+@UNPUBLISHED{Scharpf2009,

+  author = {Robert B Scharpf and Ingo Ruczinski and Benilton Carvalho and Betty

+	Doan and Aravinda Chakravarti and Rafael Irizarry},

+  title = {A multilevel model to address batch effects in copy number estimation

+	using SNP arrays},

+  month = {May},

+  year = {2009}

+}

\ No newline at end of file

@@ -13,6 +13,7 @@

 \alias{chromosome,CrlmmSetList-method}

 \alias{cnIndex,CrlmmSetList-method}

 \alias{combine,CrlmmSetList,CrlmmSetList-method}

+\alias{confs,CrlmmSetList-method}

 \alias{copyNumber,CrlmmSetList-method}

 \alias{nrow,CrlmmSetList-method}

 \alias{ncol,CrlmmSetList-method}

@@ -22,6 +23,7 @@

 \alias{sampleNames,CrlmmSetList-method}

 \alias{show,CrlmmSetList-method}

 \alias{snpIndex,CrlmmSetList-method}

+\alias{update,CrlmmSetList-method}

 \title{Class "CrlmmSetList"}

@@ -68,15 +70,15 @@

 \section{Methods}{

   \describe{

     \item{"["}{\code{signature(x = "CrlmmSetList")}: subsets both the

-  features and samples of each element in th elist}

-

+      features and samples of each element in th elist}

+  

     \item{"$"}{\code{signature(x = "CrlmmSetList")}: used to access

     element from the \code{featureData} of the \code{CopyNumberSet}

     element.  In particular, useful for extract the SNP- and

     batch-specific parameters used to estimate copy number that are

     currently stored in the \code{featureData} of \code{CopyNumberSet}

     objects. }

-  

+

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

       quantile normalized intensities for the A allele in the

       \code{ABset} element.}

@@ -84,7 +86,7 @@

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

       quantile normalized intensities for the B allele in the

       \code{ABset} element.}

-

+    

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

       batch information used to estimate copy number.}    

@@ -104,6 +106,9 @@

     \item{"combine"}{\code{signature(x = "CrlmmSetList", y = "CrlmmSetList")}: combines

       objects of the class.  }

+

+    \item{"confs"}{\code{signature(x = "CrlmmSetList")}: confidence

+    scores for crlmm genotype calls. }    

     \item{"nrow"}{\code{signature(x = "CrlmmSetList")}: Number of rows

     (features) of each element in the list.  }

@@ -125,18 +130,20 @@

     \item{"snpIndex"}{\code{signature(object = "CrlmmSetList")}: returns

       the row indices of the polymorphic probes.}

-    

-    

-  }

-}

+    \item{"update"}{\code{signature(object = "CrlmmSetList")}: This

+    method calls \code{computeCopynumber} to compute copy number for the

+    \code{crlmmSetList} object.  The object returned by this method is

+    itself an instance of class \code{CrlmmSetList}.  The third element

+    of the list is an instance of \code{CopyNumberSet} containing the

+    locus-level, allele-specific estimates of copynumber.  }

+}

+}

 \section{Warnings}{

-

    \code{combine} will produce a warning as the \code{featureData} are

    not easily combined and therefore removed.  This method needs

    improvement.

-

 }

 \author{R. Scharpf}

 \seealso{

@@ -8,7 +8,7 @@

 }

 \usage{

-computeCopynumber(object, CHR, bias.adj=FALSE, batch, SNRmin=5, cdfName="genomewidesnp6", ...)

+computeCopynumber(object, CHR, bias.adj=FALSE, batch, SNRmin=5, cdfName, ...)

 }

 \arguments{

   \item{object}{object of class \code{CrlmmSetList}.}

@@ -25,6 +25,10 @@ computeCopynumber(object, CHR, bias.adj=FALSE, batch, SNRmin=5, cdfName="genomew

 \details{

+  One may altenatively use the \code{update} method (a wrapper for

+  \code{computeCopynumber}) to add a \code{CopyNumberSet} object to the

+  \code{CrlmmSetList} object.

+