@@ -199,3 +199,24 @@ is decoded and scanned

 * fixed bugs in biasAdj and biasAdjNP

+2009-07-09 R Scharpf - committed version 1.3.10

+

+* changed defaults for crlmmWrapper: save.it=FALSE

+

+* splitByChr method looks for 'chr' or 'chromosome' in colnames 

+

+* sample-specific standard deviations in the .getEmission function

+

+* add check for duplicated positions in creating locusset object

+

+* computeCopynumber returns an object of class CrlmmSetList

+  

+  - dimnames for the elements in the list are the same

+

+  - elements in the list are ordered by chromosome and physical

+    position

+

+  - copy number parameters in the featureData of the CopyNumberSet

+    element are thresholded

+

+

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

+Version: 1.3.10

 Date: 2009-06-17

 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>

@@ -8,12 +8,12 @@ importClassesFrom(Biobase, AnnotatedDataFrame, AssayData, eSet,

 importClassesFrom(oligoClasses, SnpLevelSet)

 importMethodsFrom(Biobase, annotation, "annotation<-", annotatedDataFrameFrom,

-		  assayData, "assayData<-", combine,

+		  assayData, "assayData<-", combine, dims,

 		  experimentData, "experimentData<-",

                   fData, featureData, "featureData<-", featureNames,

 		  fvarMetadata, fvarLabels,

                   pData, phenoData, "phenoData<-", pubMedIds, rowMedians, sampleNames, scanDates,

-                  "scanDates<-", storageMode, "storageMode<-", updateObject)

+                  "scanDates<-", storageMode, "storageMode<-", updateObject, varLabels)

 ##importMethodsFrom(oligoClasses, chromosome, copyNumber, position)

@@ -24,7 +24,7 @@ importFrom(oligoClasses, chromosome, copyNumber, position, calls)

 ##importFrom(methods, "@<-", callNextMethod, new, validObject)

 importFrom(Biobase, assayDataElement, assayDataElementNames,

-           assayDataElementReplace, assayDataNew, classVersion)

+           assayDataElementReplace, assayDataNew, classVersion, validMsg)

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

            polygon, rect, segments, text, points)

@@ -46,10 +46,11 @@ importFrom(mvtnorm, dmvnorm)

 importFrom(ellipse, ellipse)

 ##S3method(ellipse,CopyNumberSet)

-

+##S3method(boxplot,CrlmmSetList)

 exportClasses(ABset, CopyNumberSet, CrlmmSetList)

 ##S3method(ellipse, CopyNumberSet)

-exportMethods(A, B,

+exportMethods("[", ##"[[",

+	      "$", A, B,

 	      calls,

 	      CA,

 	      "CA<-",

@@ -62,7 +63,8 @@ exportMethods(A, B,

 	      show,

 	      snpIndex)

-export(celDates,

+export(batch,

+       celDates,

        calls,

        chromosome,

        cnrma,

@@ -72,7 +74,7 @@ export(celDates,

        crlmm,

        crlmmIllumina,

        crlmmWrapper,

-       ellipse,

+##       ellipse,

 ##       computeEmission,

        list.celfiles,

        position,

@@ -1,9 +1,10 @@

 setGeneric("A", function(object) standardGeneric("A"))

 setGeneric("B", function(object) standardGeneric("B"))

+setGeneric("batch", function(object) standardGeneric("batch"))

 ##setGeneric("calls", function(x) standardGeneric("calls"))

 setGeneric("confs", function(object) standardGeneric("confs"))

-setGeneric("CA", function(object) standardGeneric("CA"))

-setGeneric("CB", function(object) standardGeneric("CB"))

+setGeneric("CA", function(object, ...) standardGeneric("CA"))

+setGeneric("CB", function(object, ...) standardGeneric("CB"))

 setGeneric("CA<-", function(object, value) standardGeneric("CA<-"))

 setGeneric("CB<-", function(object, value) standardGeneric("CB<-"))

 ##setGeneric("chromosome", function(object) standardGeneric("chromosome"))

@@ -11,6 +12,7 @@ setGeneric("cnIndex", function(object) standardGeneric("cnIndex"))

 setGeneric("cnNames", function(object) standardGeneric("cnNames"))

 ##setGeneric("copyNumber", function(object) standardGeneric("copyNumber"))

 ##setGeneric("position", function(object) standardGeneric("position"))

+setGeneric(".harmonizeDimnames", function(object) standardGeneric(".harmonizeDimnames"))

 setGeneric("snpIndex", function(object) standardGeneric("snpIndex"))

 setGeneric("snpNames", function(object) standardGeneric("snpNames"))

 setGeneric("splitByChromosome", function(object, ...) standardGeneric("splitByChromosome"))

@@ -206,6 +206,8 @@ harmonizeSnpSet <- function(crlmmResult, ABset){

 }

 harmonizeDimnamesTo <- function(object1, object2){

+	#object2 should be a subset of object 1

+	object2 <- object2[featureNames(object2) %in% featureNames(object1), ]

 	object1 <- object1[match(featureNames(object2), featureNames(object1)), ]

 	object1 <- object1[, match(sampleNames(object2), sampleNames(object1))]

 	stopifnot(all.equal(featureNames(object1), featureNames(object2)))

@@ -251,7 +253,8 @@ crlmmIlluminaWrapper <- function(sampleSheet, outdir="./", cdfName,

-crlmmWrapper <- function(filenames, outdir="./", cdfName="genomewidesnp6", save.it,

+crlmmWrapper <- function(filenames, outdir="./", cdfName="genomewidesnp6",

+			 save.it=FALSE,

 			 splitByChr=TRUE, ...){

 	##no visible binding for res

 	if(!file.exists(file.path(outdir, "crlmmResult.rda"))){

@@ -261,9 +264,6 @@ crlmmWrapper <- function(filenames, outdir="./", cdfName="genomewidesnp6", save.

 		message("Loading crlmmResult...")

 		load(file.path(outdir, "crlmmResult.rda"))

 	}

-

-	##- Make crlmmResult the same dimension as ABset

-	## -Create a list object that where rows and columns can be subset using '[' methods

 	if(!file.exists(file.path(outdir, "cnrmaResult.rda"))){

 		message("Quantile normalizing the copy number probes")		

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

@@ -272,23 +272,25 @@ crlmmWrapper <- function(filenames, outdir="./", cdfName="genomewidesnp6", save.

 		message("Loading cnrmaResult...")		

 		load(file.path(outdir, "cnrmaResult.rda"))

 	}

-##	loader("intensities.rda", pkgname=pkgname, envir=.crlmmPkgEnv)

-##	res <- get("res", envir=.crlmmPkgEnv)

 	load(file.path(outdir, "intensities.rda"))

 	ABset <- combineIntensities(res, cnrmaResult, cdfName)

 	scanDates(ABset) <- as.character(celDates(filenames))	

 	crlmmResult <- harmonizeSnpSet(crlmmResult, ABset)

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

-	crlmmList <- list(ABset,

-			  crlmmResult)

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

+	crlmmResults <- list(ABset,

+			     crlmmResult)

+	crlmmResults <- as(crlmmResults, "CrlmmSetList")

 	if(splitByChr){

 		message("Saving by chromosome")

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

+		splitByChromosome(crlmmResults, cdfName=cdfName, outdir=outdir)

 	} else{

 		message("Saving crlmmList object to ", outdir)

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

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

+	}

+	if(!save.it){

+		message("Cleaning up")

+		unlink(file.path(outdir, "intensities.rda"))

 	}

 	return()

 }

@@ -509,7 +511,18 @@ computeCopynumber <- function(object,

 	}

 	message("Organizing results...")			

 	locusSet <- list2locusSet(envir, ABset=ABset, NPset=NPset, CHR=CHR, cdfName=cdfName)

-	return(locusSet)

+	if(anyDuplicated(position(locusSet))){

+		message("duplicated physical positions removed from CopyNumberSet object")

+		locusSet <- locusSet[!duplicated(position(locusSet)), ]

+	}

+	message("Thresholding model parameters")

+	locusSet <- thresholdModelParams(locusSet)

+	object[[3]] <- locusSet

+	message("harmonizing dimnames of the elements in CrlmmSetList...")

+	object <- .harmonizeDimnames(object)

+	message("Reordering features by chromosome and physical position...")

+	object <- object[order(chromosome(object), position(object)), ]	

+	return(object)

 }

 cnIllumina <- function(object,

@@ -1841,21 +1854,30 @@ getParams <- function(object, batch){

 	return(params)

 }

-computeEmission <- function(crlmmResults, cnset, copyNumberStates=0:5, MIN=2^3){

+computeEmission <- function(crlmmResults, copyNumberStates=0:5, MIN=2^3,

+			    EMIT.THR,

+			    scaleSds=TRUE){

 	##threshold small nu's and phis

-	cnset <- thresholdModelParams(cnset, MIN=MIN)

+	cnset <- thresholdModelParams(crlmmResults[[3]], MIN=MIN)

 	index <- order(chromosome(cnset), position(cnset))

+	

 	if(any(diff(index) > 1)) stop("must be ordered by chromosome and physical position")

 	emissionProbs <- array(NA, dim=c(nrow(cnset), ncol(cnset), length(copyNumberStates)))

 	dimnames(emissionProbs) <- list(featureNames(crlmmResults),

 					sampleNames(crlmmResults),

 					paste("copy.number_", copyNumberStates, sep=""))	

-	batch <- cnset$batch

-	for(i in seq(along=unique(batch))){

+	b <- batch(cnset)

+	for(i in seq(along=unique(b))){

 		if(i == 1) cat("Computing emission probabilities \n")

-		message("batch ", unique(batch)[i], "...")

-		emissionProbs[, batch == unique(batch)[i], ] <- .getEmission(crlmmResults, cnset, batch=unique(batch)[i],

-				copyNumberStates=copyNumberStates)

+		message("batch ", unique(b)[i], "...")

+		emissionProbs[, b == unique(b)[i], ] <- .getEmission(crlmmResults, cnset, batch=unique(b)[i],

+				copyNumberStates=copyNumberStates,

+				scaleSds=scaleSds)

+	}

+	if(missing(EMIT.THR)){

+		EMIT.THR <- quantile(emissionProbs, probs=0.25, na.rm=TRUE)

+		message("Thresholding emission probabilities at a small negative value (", round(EMIT.THR, 1), ") to reduce influence of outliers.")

+		emissionProbs[emissionProbs < EMIT.THR] <- EMIT.THR

 	}

 	emissionProbs

 }

@@ -1877,11 +1899,35 @@ thresholdModelParams <- function(object, MIN=2^3){

 	object

 }

-.getEmission <- function(crlmmResults, cnset, batch, copyNumberStates){

+.getEmission <- function(crlmmResults, cnset, batch, copyNumberStates, scaleSds=TRUE){

 	if(length(batch) > 1) stop("batch variable not unique")

 	crlmmResults <- crlmmResults[, cnset$batch==batch]

-	cnset <- cnset[, cnset$batch == batch]	

-	emissionProbs <- array(NA, dim=c(nrow(crlmmResults[[1]]), ncol(crlmmResults[[1]]), length(copyNumberStates)))

+	cnset <- cnset[, cnset$batch == batch]

+

+##	a <- A(crlmmResults)

+##	b <- B(crlmmResults)	

+##	sds.a <- apply(log2(a), 2, sd, na.rm=TRUE)

+##	sds.b <- apply(log2(b), 2, sd, na.rm=TRUE)

+	if(scaleSds){

+		a <- CA(crlmmResults)

+		b <- CB(crlmmResults)

+		sds.a <- apply(a, 2, sd, na.rm=TRUE)

+		sds.b <- apply(b, 2, sd, na.rm=TRUE)	

+	

+		sds.a <- log2(sds.a/median(sds.a))

+		sds.b <- log2(sds.b/median(sds.b))

+		

+		sds.a <- matrix(sds.a, nrow(cnset), ncol(cnset), byrow=TRUE)

+		sds.b <- matrix(sds.b, nrow(cnset), ncol(cnset), byrow=TRUE)

+

+	} else sds.a <- sds.b <- matrix(0, nrow(cnset), ncol(cnset))

+##	ca <- CA(cnset)

+##	sds.ca <- apply(ca, 2, sd, na.rm=T)

+##	sds.ca <- sds.ca/median(sds.ca)

+##	sds.scale <- sds/median(sds)  #scale snp-specific variance by measure of the relative sample noise

+	

+	emissionProbs <- array(NA, dim=c(nrow(crlmmResults[[1]]),

+				   ncol(crlmmResults[[1]]), length(copyNumberStates)))

 	snpset <- cnset[snpIndex(cnset), ]

 	params <- getParams(snpset, batch=batch)

 	##attach(params)

@@ -1911,7 +1957,12 @@ thresholdModelParams <- function(object, MIN=2^3){

 			if(CA > 0 & CB > 0) r <- corr

 			if(CA == 0 & CB == 0) r <- 0

 			muA <- log2(nuA+CA*phiA)

-			muB <- log2(nuB+CB*phiB)		

+			muB <- log2(nuB+CB*phiB)

+

+			sigmaA <- matrix(sigmaA, nrow=length(sigmaA), ncol=ncol(cnset), byrow=FALSE)

+			sigmaB <- matrix(sigmaB, nrow=length(sigmaB), ncol=ncol(cnset), byrow=FALSE)

+			sigmaA <- sigmaA+sds.a[snpIndex(crlmmResults), ]

+			sigmaB <- sigmaB+sds.b[snpIndex(crlmmResults), ]			

 			##might want to allow the variance to be sample-specific

 			##TODO:

@@ -1954,9 +2005,10 @@ thresholdModelParams <- function(object, MIN=2^3){

 		mus <- as.numeric(matrix(log2(nuA + CT*phiA), nrow(cnset), ncol(cnset)))

 		##Again, should make sds sample-specific

 		sds.matrix <- matrix(sqrt(sig2A), nrow(cnset), ncol(cnset))

-		sds <- as.numeric(matrix(sqrt(sig2A), nrow(cnset), ncol(cnset)))

-		

-		tmp <- matrix(dnorm(a, mean=mus, sd=sds), nrow(cnset), ncol(cnset))

+

+		sds.matrix <- sds.matrix + sds.a[cnIndex(crlmmResults), ]

+		sds <- as.numeric(sds.matrix)

+		suppressWarnings(tmp <- matrix(dnorm(a, mean=mus, sd=sds), nrow(cnset), ncol(cnset)))

 		emissionProbs[cnIndex(crlmmResults), , k] <- log(tmp)

 	}

 	emissionProbs

@@ -3,8 +3,9 @@ setValidity("CopyNumberSet", function(object) {

 	msg <- validMsg(NULL, assayDataValidMembers(assayData(object), c("CA", "CB")))

 	if (is.null(msg)) TRUE else msg

 })

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

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

+##may want to allow thresholding here (... arg)

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

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

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

 		 function(object, value) assayDataElementReplace(object, "CB", value))

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

@@ -12,6 +13,14 @@ setReplaceMethod("CA", signature(object="CopyNumberSet", value="matrix"),

 setMethod("chromosome", "CopyNumberSet", function(object) fData(object)$chromosome)

+setMethod("batch", "CopyNumberSet", function(object){

+	if("batch" %in% varLabels(object)){

+		result <- object$batch

+	} else {

+		stop("batch not in varLabels of CopyNumberSet")

+	}

+	return(result)

+})

 setMethod("copyNumber", "CopyNumberSet", function(object){

 	##ensure that 2 + NA = 2 by replacing NA's with zero

@@ -29,8 +38,8 @@ setMethod("copyNumber", "CopyNumberSet", function(object){

 setMethod("position", "CopyNumberSet", function(object) fData(object)$position)

 ##setMethod("ellipse", "CopyNumberSet", function(x, copynumber, ...){

-##ellipse.CopyNumberSet <- function(x, copynumber, ...){

-setMethod("ellipse", "CopyNumberSet", function(x, copynumber, ...){

+ellipse.CopyNumberSet <- function(x, copynumber, ...){

+##setMethod("ellipse", "CopyNumberSet", function(x, copynumber, ...){

 	##fittedOrder <- unique(sapply(basename(celFiles), function(x) strsplit(x, "_")[[1]][2]))

 	##index <- match(plates, fittedOrder)

 	if(nrow(x) > 1) stop("only 1 snp at a time")

@@ -61,7 +70,7 @@ setMethod("ellipse", "CopyNumberSet", function(x, copynumber, ...){

 				      scale=scale), ...)

 		}

 	}

-})

+}

@@ -23,11 +23,28 @@ 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("A", "CrlmmSetList", function(object) A(object[[1]]))

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

+

+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("cnIndex", "CrlmmSetList", function(object) match(cnNames(object[[1]]), 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]]

@@ -52,8 +69,12 @@ setMethod("combine", signature=signature(x="CrlmmSetList", y="CrlmmSetList"),

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

+

 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"),

 	  function(x, y, ...){

 		  A <- log2(A(x))

@@ -71,18 +92,51 @@ setMethod("points", signature(x="CrlmmSetList"),

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

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

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

-	show(object[[1]])

-	show(object[[2]])

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

+})

+

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

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

+

+

+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]]), featureNames(object)))

 setMethod("splitByChromosome", "CrlmmSetList", function(object, cdfName, outdir){

 	path <- system.file("extdata", package=paste(cdfName, "Crlmm", sep=""))	

 	load(file.path(path, "snpProbes.rda"))

-	load(file.path(path, "cnProbes.rda"))				

+	load(file.path(path, "cnProbes.rda"))

+	k <- grep("chr", colnames(snpProbes))

+	if(length(k) < 1) stop("chr or chromosome not in colnames(snpProbes)")

 	for(CHR in 1:24){

 		cat("Chromosome ", CHR, "\n")

-		snps <- rownames(snpProbes)[snpProbes$chr == CHR]

-		cnps <- rownames(cnProbes)[cnProbes$chr == CHR]

+		snps <- rownames(snpProbes)[snpProbes[, k] == CHR]

+		cnps <- rownames(cnProbes)[cnProbes[, k] == CHR]

 		index <- c(match(snps, featureNames(object)),

 			   match(cnps, featureNames(object)))

 		crlmmResults <- object[index, ]

@@ -35,5 +35,27 @@ B Carvalho

     - Need accessors

+    - how to store parameters for chromosome X (should be nu_A_male,

+      nu_A_female, etc.)?

+

+  o Provide a SNP- and sample-specific estimate of the variance for

+    computing emission probabilities.  (Currently, only SNP-specific)

+

+    (DONE -- version 1.3.10

+

+  o Allow a copy-neutral ROH state by not having an equal prior on the

+    ellipses.  For instance,

+

+    normal state:  1/4 AA, 1/2 AB, 1/4 BB

+    

+    copy-neutral ROH:  (1-epsilon)/2 AA,  epsilon AB,  (1-epsilon)/2

+    BB

+

+  o add caConfidence and cbConfidence slots (but what about the correlation)

+

+

+

+

+

@@ -47,6 +47,7 @@ scanDates(example.cnset)[1:5]

 <<requiredPackages>>=

 library(crlmm)

 library(genomewidesnp6Crlmm)

+library(Biobase)

 @ 

 Specify the complete path for the CEL files and a directory in which to

@@ -111,9 +112,9 @@ batch. Here we make a table of date versus ancestry:

 <<specifyBatch>>=

 sns <- sampleNames(crlmmResults)

 sns[1]

-batch <- substr(basename(sns), 13, 13)

-table(batch)

-table(format(as.POSIXlt(scanDates(crlmmResults)), "%d %b %Y"), batch)

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

+table(plate)

+table(format(as.POSIXlt(scanDates(crlmmResults)), "%d %b %Y"), plate)

 @ 

 As all of these samples were run on the first week of March, we would

@@ -124,22 +125,35 @@ chemistry plate as an argument for \Robject{batch} in the

 \Robject{computeCopynumber} function.

 <<computeCopynumber>>=

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

-	if(file.exists(file.path(outdir, paste("cnset_", CHR, ".rda", sep="")))) load(file.path(outdir, paste("cnset_", CHR, ".rda", sep="")))	

-	else {

-		cnset <- computeCopynumber(crlmmResults, SNRmin=5, batch=batch, CHR=CHR, cdfName="genomewidesnp6")

-		save(cnset, file=file.path(outdir, paste("cnset_", CHR, ".rda", sep="")))

-	}

+if(!exists("crlmmResults2")){

+	crlmmResults2 <- crlmmResults

+	crlmmResults <- computeCopynumber(crlmmResults, SNRmin=5, batch=plate, CHR=CHR, cdfName="genomewidesnp6")

+	scanDates(crlmmResults[[3]]) <- scanDates(crlmmResults[[3]])

 }

-scanDates(cnset) <- scanDates(crlmmResults)

-cnset <- cnset[order(chromosome(cnset), position(cnset)), ]

-crlmmResults <- crlmm:::harmonizeDimnamesTo(crlmmResults, cnset)

 @ 

+**Note: the number of samples in the \Robject{CrlmmSetList} object after

+copy number estimation may be fewer than the number of samples in the

+\Robject{CrlmmSetList} object after preprocessing/genotyping.  This

+occurs when 1 or more samples have a signal-to-noise ratio less than

+value passed to \Robject{SNRmin}.  By default, intermediate forms of the

+data are stored in one object to ensure that each element in the

+\Robject{CrlmmSetList} have the same ordering of probes and samples.

+The object returned by \Robject{computeCopynumber} is ordered by

+chromosome and physical position (useful for downstream methods that

+smooth the copy number as a function of the physical position).  **

+

+<<dimnamesMayDiffer>>=

+identical(featureNames(crlmmResults), featureNames(crlmmResults2)) ##not identical

+identical(sampleNames(crlmmResults), sampleNames(crlmmResults2))

+gc()

+@ 

+

+

 <<example.cnset, eval=FALSE, echo=FALSE>>=

 example.crlmmResults <- crlmmResults[1:1000, 1:100]

-example.cnset <- cnset[1:1000, 1:100]

-save(example.cnset, file="~/madman/Rpacks/crlmm/inst/scripts/example.cnset.rda")

+##example.cnset <- cnset[1:1000, 1:100]

+##save(example.cnset, file="~/madman/Rpacks/crlmm/inst/scripts/example.cnset.rda")

 save(example.crlmmResults, file="~/madman/Rpacks/crlmm/inst/scripts/example.crlmmResults.rda")

 @ 

@@ -151,8 +165,8 @@ additional robustness to this assumption.  Set the \Robject{bias.adj}

 argument to \Robject{TRUE}:

 <<biasAdjustment, eval=FALSE>>=

-cnset <- computeCopynumber(crlmmResults, SNRmin=5, batch=batch, CHR=CHR, 

-			   cdfName="genomewidesnp6", bias.adj=TRUE)

+crlmmResults[[3]] <- computeCopynumber(crlmmResults, SNRmin=5, batch=plate, CHR=CHR, 

+				       cdfName="genomewidesnp6", bias.adj=TRUE)

 @ 

 \section{Suggested plots}

@@ -165,12 +179,12 @@ integer.

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

 par(las=1)

-plot(position(cnset), copyNumber(cnset)[, 1], pch=".", cex=2, xaxt="n", col="grey20", ylim=c(0,6), 

+plot(position(crlmmResults), copyNumber(crlmmResults)[, 1], pch=".", cex=2, xaxt="n", col="grey20", ylim=c(0,6), 

      ylab="copy number", xlab="physical position (Mb)",

-     main=paste(sampleNames(cnset)[1], ", CHR:", unique(chromosome(cnset))))

-points(position(cnset)[cnIndex(cnset)], copyNumber(cnset)[cnIndex(cnset), 1],

+     main=paste(sampleNames(crlmmResults)[1], ", CHR:", unique(chromosome(crlmmResults))))

+points(position(crlmmResults)[cnIndex(crlmmResults)], copyNumber(crlmmResults)[cnIndex(crlmmResults), 1],

        pch=".", cex=2, col="lightblue")

-axis(1, at=pretty(range(position(cnset))), labels=pretty(range(position(cnset)))/1e6)

+axis(1, at=pretty(range(position(crlmmResults))), labels=pretty(range(position(crlmmResults)))/1e6)

 @ 

 <<idiogram, eval=FALSE, echo=FALSE>>=

@@ -216,7 +230,7 @@ for(i in indices[[j]]){

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

 require(RColorBrewer)

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

-J <- split(1:ncol(cnset), cnset$batch)

+J <- split(1:ncol(crlmmResults), batch(crlmmResults))

 colors <- c("red", "blue", "green3")

 cex <- 0.6

 colors <- c("blue", greens[8], "red")

@@ -242,15 +256,15 @@ for(j in seq(along=indices)[1:10]){

 		     xlim=xlim, ylim=ylim,

 		     type="n")

 		if(plotpoints){

-			for(b in seq(along=unique(cnset$batch))){

+			for(b in seq(along=unique(batch(crlmmResults)))){

 				points(crlmmResults[i, J[[b]]], 

 				       pch=pch, 

 				       col=colors[b], bg=colors[b], cex=cex,

 				       xlim=xlim, ylim=ylim)

 			}

 		}

-		for(b in seq(along=unique(cnset$batch))){

-			ellipse(cnset[i, J[[b]]], copynumber=2, col=colors[b], lwd=lwd)

+		for(b in seq(along=unique(batch(crlmmResults)))){

+			ellipse(crlmmResults[i, J[[b]]], copynumber=2, col=colors[b], lwd=lwd)

 		}

 		##legend("bottomright", bty="n", legend=featureNames(crlmmResults)[i])

 		if(k == 1) {

@@ -264,10 +278,7 @@ for(j in seq(along=indices)[1:10]){

 dev.off()

 @

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

-##8300117, 8469401, 8706404, 2296864, 8659838, 1936806, 4198761, 8584395

-@ 

 \begin{figure}

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

@@ -285,14 +296,13 @@ with the Birdseye HMM, the same transition probabilities are specified.

 <<segment>>=

 library(VanillaICE)

 copyNumberStates <- 0:5

-require(Biobase)

-cnset <- crlmm:::thresholdModelParams(cnset)

+##crlmmResults <- crlmm:::thresholdModelParams(crlmmResults)

 if(!exists("fit")){

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

 	else {

-		emission.cn <- computeEmission(crlmmResults, cnset, copyNumberStates)

-		tau <- transitionProbability(chromosome=chromosome(cnset),

-					     position=position(cnset),

+		emission.cn <- crlmm:::computeEmission(crlmmResults, copyNumberStates)

+		tau <- transitionProbability(chromosome=chromosome(crlmmResults),

+					     position=position(crlmmResults),

 					     TAUP=1e8)

 		initialP <- rep(1/length(copyNumberStates), length(copyNumberStates))

 		fit <- viterbi(initialStateProbs=log(initialP),

@@ -304,8 +314,7 @@ if(!exists("fit")){

 			       altered2normal=0.5,

 			       altered2altered=0.0025)

 		brks <- breaks(x=fit, states=copyNumberStates, position=tau[, "position"],

-			       chromosome=tau[, "chromosome"],

-			       sampleNames=sampleNames(cnset))

+			       chromosome=tau[, "chromosome"])

 	}

 }

 @ 

@@ -314,31 +323,27 @@ if(!exists("fit")){

 dir.create("figures")

 par(las=1, mfrow=c(1,1), ask=FALSE)

 bmp("figures/chr22plots%02d.png", width=800, height=500)

-for(j in 1:ncol(cnset)){

+for(j in 1:ncol(crlmmResults)){

 	cat(j, " ")

 	##A region that we would miss if we only looked at the polymorphic probes...is it real?

-	##index <- which(position(cnset) > 21.2*1e6 & position(cnset) < 21.6*1e6 & mns < 1)

+	##index <- which(position(crlmmResults) > 21.2*1e6 & position(crlmmResults) < 21.6*1e6 & mns < 1)

 	##fit <- fit[index,]

 	index <- 1:nrow(fit)

-	plot(0:1, type="n", xlim=range(position(cnset)[index]),

+	plot(0:1, type="n", xlim=range(position(crlmmResults)[index]),

 	     xaxt="n", ylim=c(-0.5, 6), 

 	     ylab="copy number", xlab="physical position (Mb)",

-	     main=paste(sampleNames(cnset)[j], ", CHR:", unique(chromosome(cnset))), 

+	     main=paste(sampleNames(crlmmResults)[j], ", CHR:", unique(chromosome(crlmmResults))), 

 	     lwd=2,

 	     col="blue")

-	points(position(cnset)[index], copyNumber(cnset)[index, j], pch=".", cex=2, col="grey60")

-	lines(position(cnset)[index], fit[index, j]-1, type="s", lwd=2, col="blue")

+	points(position(crlmmResults)[index], copyNumber(crlmmResults)[index, j], pch=".", cex=2, col="grey60")

+	lines(position(crlmmResults)[index], fit[index, j]-1, type="s", lwd=2, col="blue")

 	abline(h=0:4, col="grey60")

-	axis(1, at=pretty(range(position(cnset)[index])), labels=pretty(range(position(cnset)[index]))/1e6)

+	axis(1, at=pretty(range(position(crlmmResults)[index])), labels=pretty(range(position(crlmmResults)[index]))/1e6)

 }

 dev.off()

 @ 

-<<fixProblems, eval=FALSE>>=

-## NA18522, NA18978 (missed deletion), NA18944 (noisy), NA19098 and NA19130 (missed amplification), 

-## NA19138 (looks like deletion with contamination), NA19143 (missed amplification)

-## NA19204 (looks like homozygous deletion)

-@ 

+

 \section{Session information}

 <<sessionInfo, results=tex>>=

@@ -346,4 +351,6 @@ toLatex(sessionInfo())

 @ 

+

+

 \end{document}

@@ -2,6 +2,7 @@

 \Rdversion{1.1}

 \docType{class}

 \alias{CopyNumberSet-class}

+\alias{batch,CopyNumberSet-method}

 \alias{CA}

 \alias{CA<-}

 \alias{CA,CopyNumberSet-method}

@@ -40,14 +41,17 @@ Class \code{"\linkS4class{Versioned}"}, by class "eSet", distance 3.

 \section{Methods}{

   \describe{

-    \item{CA}{\code{signature(object = "CopyNumberSet")}: Extracts the

+    \item{batch}{\code{signature(object="CopyNumberSet")}: Extracts the

+    batch information used to estimate copy number.}

+

+    \item{CA}{\code{signature(object = "CopyNumberSet", ...)}: Extracts the

       copy number for allele 'A' at polymorphic loci. At nonpolymorphic

       loci, the total copy number is returned.}

     \item{CA<-}{\code{signature(object = "CopyNumberSet",

     value="matrix")}: Replaces CA matrix with supplied value.}    

-    \item{CB}{\code{signature(object = "CopyNumberSet")}:

+    \item{CB}{\code{signature(object = "CopyNumberSet", ...)}:

       Extracts copy number for allele 'B' at polymorphic loci. NAs are

       returned for nonpolymorphic loci.}    

@@ -3,13 +3,22 @@

 \docType{class}

 \alias{CrlmmSetList-class}

 \alias{[,CrlmmSetList-method}

+\alias{$,CrlmmSetList-method}

 \alias{A,CrlmmSetList-method}

 \alias{B,CrlmmSetList-method}

+\alias{batch,CrlmmSetList-method}

+\alias{CA,CrlmmSetList-method}

+\alias{CB,CrlmmSetList-method}

 \alias{calls,CrlmmSetList-method}

+\alias{chromosome,CrlmmSetList-method}

 \alias{cnIndex,CrlmmSetList-method}

 \alias{combine,CrlmmSetList,CrlmmSetList-method}

+\alias{copyNumber,CrlmmSetList-method}

+\alias{nrow,CrlmmSetList-method}

+\alias{ncol,CrlmmSetList-method}

 \alias{plot,CrlmmSetList,ANY-method}

 \alias{points,CrlmmSetList-method}

+\alias{position,CrlmmSetList-method}

 \alias{sampleNames,CrlmmSetList-method}

 \alias{scanDates,CrlmmSetList-method}

 \alias{show,CrlmmSetList-method}

@@ -59,17 +68,35 @@

 }

 \section{Methods}{

   \describe{

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

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

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

-    

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

+

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

+  

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

+      copy number for allele B at polymorphic loci. For nonpolymorphic

+      probes, CB returns 'NA'.}    

+

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

       the genotype calls from the \code{SnpSet} element.}

@@ -78,7 +105,13 @@

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

       objects of the class.  }

-

+    

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

+    (features) of each element in the list.  }

+    

+    \item{"ncol"}{\code{signature(x = "CrlmmSetList")}: Number of columns

+    (samples) of each element in the list.  }    

+	

     \item{"plot"}{\code{signature(x = "CrlmmSetList")}: For A versus B

       scatterplots.  }

new file mode 100644

@@ -0,0 +1,25 @@

+\name{batch}

+\Rdversion{1.1}

+\alias{batch}

+\title{

+  Stores information on the batch (e.g., chemistry plate) used to

+  process samples.

+}

+\description{

+  Parameters for estimating copy number are locus- and batch-specific.

+}

+\usage{

+batch(object)

+}

+\arguments{

+  \item{object}{An object of class \code{CrlmmSetList} or \code{CopyNumberSet}.

+  }

+}

+\value{

+  Typically, a character string or factor.

+}

+\author{

+  Rob Scharpf

+}

+\keyword{manip}

+

@@ -10,7 +10,7 @@

   confidence scores are assigned using the crlmm algorithm.

 }

 \usage{

-crlmmWrapper(filenames, outdir = "./", cdfName = "genomewidesnp6", save.it, splitByChr=TRUE, ...)

+crlmmWrapper(filenames, outdir = "./", cdfName = "genomewidesnp6", save.it=FALSE, splitByChr=TRUE, ...)

 }

 \arguments{

   \item{filenames}{