@@ -174,6 +174,9 @@ is decoded and scanned

 * several changes to the copynumber vignette

+2009-06-29 R Scharpf - committed version 1.3.6

+* updated computeCopynumber

+* export A, B methods for ABset class

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

+Version: 1.3.6

 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>

@@ -45,9 +45,12 @@ importFrom(mvtnorm, dmvnorm)

 importFrom(ellipse, ellipse)

+S3method(ellipse,CopyNumberSet)

+

 exportClasses(ABset, CopyNumberSet, CrlmmSetList)

 ##S3method(ellipse, CopyNumberSet)

-exportMethods(calls,

+exportMethods(A, B,

+	      calls,

 	      CA,

 	      "CA<-",

 	      CB,

new file mode 100644

@@ -0,0 +1,18 @@

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

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

+##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, value) standardGeneric("CA<-"))

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

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

+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("snpIndex", function(object) standardGeneric("snpIndex"))

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

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

+

@@ -456,8 +456,7 @@ computeCopynumber <- function(object,

 			   bias.adj=FALSE,

 			   SNRmin=SNRmin,

 			   ...)

-	message("Organizing results...")	

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

+

 	if(bias.adj){

 		message("Running bias adjustment...")

 		.computeCopynumber(chrom=CHR,

@@ -472,10 +471,9 @@ computeCopynumber <- function(object,

 				   bias.adj=TRUE,

 				   SNRmin=SNRmin,

 				   ...)

-		message("Organizing results...")			

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

-		##.GlobalEnv[["locusSetAdj"]] <- locusSetAdj

 	}

+	message("Organizing results...")			

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

 	return(locusSet)

 }

new file mode 100644

@@ -0,0 +1,11 @@

+setValidity("ABset", function(object) {

+	##msg <- validMsg(NULL, Biobase:::isValidVersion(object, "CopyNumberSet"))

+	msg <- validMsg(NULL, assayDataValidMembers(assayData(object), c("A", "B")))

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

+})

+setMethod("A", "ABset", function(object) assayData(object)[["A"]])

+setMethod("B", "ABset", function(object) assayData(object)[["B"]])

+

+

+

+

new file mode 100644

@@ -0,0 +1,65 @@

+setValidity("CopyNumberSet", function(object) {

+	##msg <- validMsg(NULL, Biobase:::isValidVersion(object, "CopyNumberSet"))

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

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

+})

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

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

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

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

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

+		 function(object, value) assayDataElementReplace(object, "CA", value))

+

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

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

+

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

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

+	CA <- CA(object)

+	CB <- CB(object)

+	nas <- is.na(CA) & is.na(CB)

+	CA[is.na(CA)] <- 0

+	CB[is.na(CB)] <- 0

+	CN <- CA/100 + CB/100

+	##if both CA and CB are NA, report NA

+	CN[nas] <- NA

+	CN

+})

+

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

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

+	batch <- unique(x$batch)

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

+	nuA <- as.numeric(fData(x)[, match(paste("nuA", batch, sep="_"), fvarLabels(x))])

+	nuB <- as.numeric(fData(x)[, match(paste("nuB", batch, sep="_"), fvarLabels(x))])	

+	phiA <- as.numeric(fData(x)[, match(paste("phiA", batch, sep="_"), fvarLabels(x))])

+	phiB <- as.numeric(fData(x)[, match(paste("phiB", batch, sep="_"), fvarLabels(x))])	

+	tau2A <- as.numeric(fData(x)[, match(paste("tau2A", batch, sep="_"), fvarLabels(x))])

+	tau2B <- as.numeric(fData(x)[, match(paste("tau2B", batch, sep="_"), fvarLabels(x))])

+	sig2A <- as.numeric(fData(x)[, match(paste("sig2A", batch, sep="_"), fvarLabels(x))])

+	sig2B <- as.numeric(fData(x)[, match(paste("sig2B", batch, sep="_"), fvarLabels(x))])

+	corrA.BB <- as.numeric(fData(x)[, match(paste("corrA.BB", batch, sep="_"), fvarLabels(x))])

+	corrB.AA <- as.numeric(fData(x)[, match(paste("corrB.AA", batch, sep="_"), fvarLabels(x))])

+	corr <- as.numeric(fData(x)[, match(paste("corr", batch, sep="_"), fvarLabels(x))])

+	for(CN in copynumber){

+		for(CA in 0:CN){

+			CB <- CN-CA

+			A.scale <- sqrt(tau2A*(CA==0) + sig2A*(CA > 0))

+			B.scale <- sqrt(tau2B*(CB==0) + sig2B*(CB > 0))

+			scale <- c(A.scale, B.scale)

+			if(CA == 0 & CB > 0) rho <- corrA.BB

+			if(CA > 0 & CB == 0) rho <- corrB.AA

+			if(CA > 0 & CB > 0) rho <- corr

+			lines(ellipse(x=rho, centre=c(log2(nuA+CA*phiA),

+					     log2(nuB+CB*phiB)),

+				      scale=scale), ...)

+		}

+	}

+}

+

+

+

+

new file mode 100644

@@ -0,0 +1,90 @@

+setMethod("[", "CrlmmSetList", function(x, i, j, ..., drop = FALSE){

+            if (missing(drop)) drop <- FALSE

+            if (missing(i) && missing(j))

+              {

+                 if (length(list(...))!=0)

+                   stop("specify genes or samples to subset; use '",

+                        substitute(x), "$", names(list(...))[[1]],

+                        "' to access phenoData variables")

+                 return(x)

+               }

+            if (!missing(j)){

+              f1 <- function(x, j){

+                x <- x[, j]

+              }

+              x <- lapply(x, f1, j)

+            }

+            if(!missing(i)){

+              f2 <- function(x, i){

+                x <- x[i, ]

+              }

+              x <- lapply(x, f2, i)

+            }

+	as(x, "CrlmmSetList")	

+})

+

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

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

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

+setMethod("cnIndex", "CrlmmSetList", function(object) match(cnNames(object[[1]]), featureNames(object)))

+

+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("removing featureData")

+		  ##fd1 <- featureData(x.snpset)

+		  ##fd2 <- featureData(y.snpset)

+		  featureData(x.snpset) <- annotatedDataFrameFrom(calls(x.snpset), byrow=TRUE)

+		  featureData(y.snpset) <- annotatedDataFrameFrom(calls(y.snpset), byrow=TRUE)

+		  snpset <- combine(x.snpset, y.snpset)

+		  merged <- list(abset, snpset)

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

+		  merged

+	  })

+		  

+

+

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

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

+	  function(x, y, ...){

+		  A <- log2(A(x))

+		  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("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]])

+})

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

+	for(CHR in 1:24){

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

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

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

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

+			   match(cnps, featureNames(object)))

+		crlmmResults <- object[index, ]

+		save(crlmmResults, file=file.path(outdir, paste("crlmmResults_", CHR, ".rda", sep="")))

+	}

+})

+

+

new file mode 100644

@@ -0,0 +1,31 @@

+setMethod("cnIndex", "eSet", function(object) match(cnNames(object), featureNames(object)))

+setMethod("cnNames", "eSet", function(object) featureNames(object)[grep("CN_", featureNames(object))])

+##setMethod("combine", signature=signature(x="eSet", y="eSet"),

+##	  function(x, y, ...){

+##		  ##Check that both x and y are valid objects

+##		  if(!validObject(x)) stop("x is not a valid object")

+##		  if(!validObject(y)) stop("y is not a valid object")

+##

+##		  if(annotation(x) != annotation(y)) stop("must have same value for annotation slot")

+##		  if(class(x) != class(y)) stop("objects must have the same class")

+##		  if(storageMode(assayData(x)) != storageMode(assayData(y))){

+##			  stop("objects must have same storage mode for assayData")

+##		  }

+##		  fd <- combine(featureData(x), featureData(y))

+##		  pd <- combine(phenoData(x), phenoData(y))            

+##		  ad.x <- as.list(assayData(x))

+##		  ad.y <- as.list(assayData(y))

+##		  ad.xy <- mapply(rbind, ad.x, ad.y, SIMPLIFY=FALSE)

+##		  id.x <- match(rownames(ad.xy[[1]]), featureNames(fd))

+##		  ee <- combine(experimentData(x), experimentData(y))

+##		  assayData(x) <- ad.xy

+##		  storageMode(assayData(x)) <- storageMode(assayData(y))            

+##		  experimentData(x) <- ee

+##		  featureData(x) <- fd

+##		  phenoData(x) <- pd

+##		  x

+##          })

+setMethod("snpIndex", "eSet", function(object) match(snpNames(object), featureNames(object)))

+setMethod("snpNames", "eSet", function(object) featureNames(object)[grep("SNP_", featureNames(object))])

+

+

@@ -19,16 +19,21 @@ B Carvalho

 - crlmm should return results in CHP files

-#####################################

-### FOR CNRMA

-#####################################

-- Bias adjustment for X chromosome

+  ***********************************************************

+  *                                                         *       

+  *                FOR CNRMA                                *

+  *                                                         *

+  ***********************************************************

+

+  o Bias adjustment for X chromosome

+

+  o Adjust nu for altered copy number (crosshyb)

-- adjust nu for altered copy number (crosshyb)

+  o Should store parameters in something besides the featureData slot

-- should store parameters in something besides the featureData slot

+    - Perhaps add a matrix slot for parameters and a slot for batch.

-        - Perhaps add a slot for parameters and a slot for batch.

+    - Need accessors

@@ -146,24 +146,13 @@ save(example.crlmmResults, file="~/madman/Rpacks/crlmm/inst/scripts/example.crlm

 The above algorithm for estimating copy number is predicated on the

 assumption that most samples within a batch have copy number 2 at any

 given locus.  For common copy number variants, this assumption may not

-hold.  An additional iteration using a bias correction can improve the

-estimates.  Set the \Robject{bias.adj} argument to \Robject{TRUE}.

-

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

-library(oligoClasses);library(genefilter)

-source("~/madman/Rpacks/crlmm/R/AllClasses.R")

-source("~/madman/Rpacks/crlmm/R/AllGenerics.R")

-source("~/madman/Rpacks/crlmm/R/cnrma-functions.R")

-source("~/madman/Rpacks/crlmm/R/methods-eSet.R")

-source("~/madman/Rpacks/crlmm/R/methods-ABset.R")

-source("~/madman/Rpacks/crlmm/R/methods-CrlmmSetList.R")

-source("~/madman/Rpacks/crlmm/R/methods-CopyNumberSet.R")

-source("~/madman/Rpacks/crlmm/R/utils.R")

-.crlmmPkgEnv <- new.env()

-envir2 <- envir

-envir <- getCopyNumberEnvironment(crlmmSetList, cnSet)

-updateNuPhi(crlmmSetList, cnSet)

-##TODO

+hold.  An additional iteration using a bias correction provides

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

 @ 

 \section{Suggested plots}

@@ -352,8 +341,8 @@ dev.off()

 @ 

 \section{Session information}

-<<>>=

-sessionInfo()

+<<sessionInfo, results=tex>>=

+toLatex(sessionInfo())

 @ 

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

@@ -46,7 +46,9 @@ if(!exists("res.rda")){

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

 	}

 }

-@ 

+@

+

+TODO: The code below needs to be udpated for the Illumina platform.

 <<SNR>>=

 hist(crlmmOut$SNR) ##approx. 5-fold higher than what we see in Affy!