@@ -60,7 +60,15 @@ BioC data packages

 * fixed documentation to acommodate this change

 2009-03-28 B Carvalho - committed version 1.0.64

+* modified crlmm and cnrma to look for data at extdata/

-* modified warning message

+2009-03-29 Rob Scharpf - committed version 1.0.65

-* modified crlmm and cnrma to look for data at extdata/

+* update to copynumber vignette

+

+* updated computeCopynumber.Rd to accommodate SnpSet objects and

+  changes to genomewidesnp6Crlmm

+

+* cnrma function requires the cdfName

+

+* update man files for cnrma and computeCopynumber

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

 Package: crlmm

 Type: Package

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

-Version: 1.0.64

+Version: 1.0.65

 Date: 2008-12-28

 Author: Rafael A Irizarry, Benilton S Carvalho <bcarvalh@jhsph.edu>, Robert Scharpf <rscharpf@jhsph.edu>

 Maintainer: Benilton S Carvalho <bcarvalh@jhsph.edu>, Robert Scharpf <rscharpf@jhsph.edu>

@@ -28,6 +28,7 @@ generateIXTX <- function(x, nrow=3) {

 	XTX <- crossprod(X)

 	solve(XTX)

 }

+

 nuphiAllele <- function(p, allele, Ystar, W, envir){

 	Ns <- envir[["Ns"]]

 	CHR <- envir[["chrom"]]

@@ -120,7 +121,8 @@ celDates <- function(celfiles){

 	return(celdts)

 }

-cnrma <- function(filenames, sns, cdfName, seed=1, verbose=FALSE){

+cnrma <- function(filenames, cdfName="genomewidesnp6", sns, seed=1, verbose=FALSE){

+	if(cdfName != "genomewidesnp6") stop("Only genomewidesnp6 supported at this time")

 	## BC: 03/14/09

 	## getting pkgname from cdfName, in the future this might be useful

 	## as the method might be implemented for other platforms

@@ -131,7 +133,6 @@ cnrma <- function(filenames, sns, cdfName, seed=1, verbose=FALSE){

         loader("npProbesFid.rda", .crlmmPkgEnv, pkgname)

 ##	data("npProbesFid", package=pkgname, envir=.crlmmPkgEnv)

 	fid <- getVarInEnv("npProbesFid")

-	gc()

 	set.seed(seed)

 	idx2 <- sample(length(fid), 10^5) ##for skewness. no need to do everything

 	SKW <- vector("numeric", length(filenames))

@@ -378,7 +379,7 @@ computeCopynumber <- function(chrom,

 				    envir=envir)

 		}

 		steps[3] <- TRUE

-		envir[["snp-cn"]] <- steps						

+		envir[["steps"]] <- steps						

 	}

 	if(!steps[4]){

 		message("\nCopy number for nonpolymorphic probes...")	

@@ -389,7 +390,7 @@ computeCopynumber <- function(chrom,

 				       envir=envir)

 		}

 		steps[4] <- TRUE

-		envir[["np-cn"]] <- steps

+		envir[["step"]] <- steps

 	}

 }

@@ -1053,14 +1054,99 @@ biasAdj <- function(plateIndex, envir, priorProb){

 }

+posteriorNonpolymorphic <- function(plateIndex, envir, priorProb, cnStates=0:6){

+	p <- plateIndex

+	CHR <- envir[["chrom"]]

+	if(missing(priorProb)) priorProb <- rep(1/length(cnStates), length(cnStates)) ##uniform	

+	plate <- envir[["plate"]]

+	uplate <- envir[["plate"]]

+	NP <- envir[["NP"]][, plate==uplate[p]]

+	nuT <- envir[["nuT"]][, p]

+	phiT <- envir[["phiT"]][, p]

+	sig2T <- envir[["sig2T"]][, p]

+	##Assuming background variance for np probes is the same on the log-scale

+	emit <- array(NA, dim=c(nrow(NP), ncol(NP), length(cnStates)))##SNPs x sample x 'truth'

+	lT <- log2(NP)

+	sds <- sqrt(sig2T)

+	counter <- 1##state counter	

+	for(CT in cnStates){

+		cat(".")

+		if(CHR == 23) browser()

+		means <- suppressWarnings(log2(nuT + CT*phiT))

+		emit[, , counter] <- dnorm(lT, mean=means, sd=sds)

+		counter <- counter+1

+	}

+	for(j in seq(along=cnStates)){

+		emit[, , j] <- priorProb[j]*emit[, , j]

+	}

+	homDel <- emit[, , 1]

+	hemDel <- emit[, , 2]

+	norm <- emit[, , 3]

+	amp <- emit[, , 4]

+	amp4 <- emit[, , 5]

+	amp5 <- emit[, , 6]

+	amp6 <- emit[, , 7]

+	total <- homDel+hemDel+norm+amp+amp4+amp5+amp6

+	weights <- array(NA, dim=c(nrow(NP), ncol(NP), length(cnStates)))

+	weights[, , 1] <- homDel/total

+	weights[, , 2] <- hemDel/total

+	weights[, , 3] <- norm/total

+	weights[, , 4] <- amp/total

+	weights[, , 5] <- amp4/total

+	weights[, , 6] <- amp5/total

+	weights[, , 7] <- amp6/total

+	##posterior mode

+	posteriorMode <- apply(weights, c(1, 2), function(x) order(x, decreasing=TRUE)[1])

+	posteriorMode <- posteriorMode-1

+	##sns <- envir[["sns"]]

+	##colnames(posteriorMode) <- sns

+	##envir[["np.posteriorMode"]] <- posteriorMode

+	##envir[["np.weights"]] <- weights

+	posteriorMeans <- 0*homDel/total + 1*hemDel/total + 2*norm/total + 3*amp/total + 4*amp4/total + 5*amp5/total + 6*amp6/total

+	##colnames(posteriorMeans) <- sns

+	##envir[["np.posteriorMeans"]] <- posteriorMeans

+	return(posteriorMode)

+}

+

+posteriorWrapper <- function(envir){

+	snp.PM <- matrix(NA, length(envir[["snps"]]), length(envir[["sns"]]))

+	np.PM <- matrix(NA, length(envir[["cnvs"]]), length(envir[["sns"]]))

+	plate <- envir[["plate"]]

+	uplate <- envir[["uplate"]]

+	for(p in seq(along=uplate)){

+		tmp <- expectedC(plateIndex=p, envir=envir)

+		snp.PM[, plate==uplate[p]] <- tmp

+		##snp.pm <- env[["posteriorMode"]]

+		##trace(posteriorNonpolymorphic, browser)

+		tmp <- posteriorNonpolymorphic(plateIndex=p, envir=envir)

+		np.PM[, plate==uplate[p]] <- tmp##env[["np.posteriorMode"]]

+		##pMode <- rbind(snp.pm, np.pm)

+		##rownames(pMode) <- c(env[["snps"]], env[["cnvs"]])

+		##dn <- dimnames(pMode)

+		##pMode <- matrix(as.integer(pMode), nrow(pMode), ncol(pMode))

+	}

+	PM <- rbind(snp.PM, np.PM)

+	PM <- matrix(as.integer(PM), nrow(PM), ncol(PM))

+	dns <- list(c(envir[["snps"]], envir[["cnvs"]]), envir[["sns"]])

+	dimnames(PM) <- dns

+	return(PM)

+}

+

+

+

 ##for polymorphic probes

-expectedC <- function(A, B, plateIndex, envir, priorProb){

+expectedC <- function(plateIndex, envir, priorProb, cnStates=0:6){

 	p <- plateIndex

 	CHR <- envir[["chrom"]]

-	if(missing(priorProb)) priorProb <- rep(1/4, 6) ##uniform	

+	if(missing(priorProb)) priorProb <- rep(1/length(cnStates), length(cnStates)) ##uniform	

 	plate <- envir[["plate"]]

+	uplate <- envir[["uplate"]]

+	A <- envir[["A"]]

+	B <- envir[["B"]]

+	A <- A[, plate==uplate[p]]

+	B <- B[, plate==uplate[p]]

 	calls <- envir[["calls"]]	

 	calls <- calls[, plate==unique(plate)[p]]

 	probA <- sqrt(rowMeans(calls == 1, na.rm=TRUE))

@@ -1084,7 +1170,7 @@ expectedC <- function(A, B, plateIndex, envir, priorProb){

 	lB <- log2(B)	

 	X <- cbind(lA, lB)	

 	counter <- 1##state counter

-	for(CT in 0:6){

+	for(CT in cnStates){

 		cat(".")

 		for(CA in 0:CT){

 			CB <- CT-CA

@@ -1122,9 +1208,8 @@ expectedC <- function(A, B, plateIndex, envir, priorProb){

 	norm <- priorProb[3]*emit[, , 4:6]

 	amp <- priorProb[4]*emit[, , 7:10]

 	amp4 <- priorProb[5]*emit[, , 11:15]

-	amp5 <- priorProb[5]*emit[, , 16:21]

-	amp6 <- priorProb[5]*emit[, , 22:28]	

-	

+	amp5 <- priorProb[6]*emit[, , 16:21]

+	amp6 <- priorProb[7]*emit[, , 22:28]	

 	##sum over the different combinations within each copy number state

 	hemDel <- apply(hemDel, c(1,2), sum)

 	norm <- apply(norm, c(1, 2), sum)

@@ -1144,18 +1229,24 @@ expectedC <- function(A, B, plateIndex, envir, priorProb){

 	##posterior mode

 	posteriorMode <- apply(weights, c(1, 2), function(x) order(x, decreasing=TRUE)[1])

 	posteriorMode <- posteriorMode-1

-	envir[["posteriorMode"]] <- posteriorMode

+	##This is for one plate.  Need to instantiate a much bigger

+	##object in the environment

+	

+	##envir[["posteriorMode"]] <- posteriorMode

 	##weights <- list(homDel/total, hemDel/total, norm/total, amp/total, amp4/total, amp5/total, amp6/total)

-	##names(weights) <- c(0:6)

-	envir[["weights"]] <- weights

+	##names(weights) <- c(cnStates)

+	##envir[["weights"]] <- weights

 	posteriorMeans <- 0*homDel/total + 1*hemDel/total + 2*norm/total + 3*amp/total + 4*amp4/total + 5*amp5/total + 6*amp6/total

-	sns <- envir[["sns"]]

-	colnames(posteriorMeans) <- sns

-	envir[["posteriorMeans"]] <- posteriorMeans

+	##sns <- envir[["sns"]]

+	##colnames(posteriorMeans) <- sns

+	##envir[["posteriorMeans"]] <- posteriorMeans

+	return(posteriorMode)

 }

+

+

 biasAdjNP <- function(plateIndex, envir, priorProb){

 	p <- plateIndex

 	normalNP <- envir[["normalNP"]]

@@ -136,11 +136,11 @@ list2SnpSet <- function(x, returnParams=FALSE){

 }

 loader <- function(theFile, envir, pkgname){

-  stopifnot(theFile %in% c("genotypeStuff.rda", "mixtureStuff.rda", "preprocStuff.rda"))

-  theFile <- file.path(system.file(package=pkgname),

-                       "extdata", theFile)

-  if (!file.exists(theFile))

-    stop("File", theFile, "does not exist in", pkgname)

-  load(theFile, envir=envir)

+	##stopifnot(theFile %in% c("genotypeStuff.rda", "mixtureStuff.rda", "preprocStuff.rda"))

+	theFile <- file.path(system.file(package=pkgname),

+			     "extdata", theFile)

+	if (!file.exists(theFile))

+		stop("File", theFile, "does not exist in", pkgname)

+	load(theFile, envir=envir)

 }

@@ -70,7 +70,7 @@ Quantile normalize the nonpolymorphic probes and save the results.

 if(!exists("cnrmaResult")){

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

 	else {

-		cnrmaResult <- cnrma(celFiles)

+		cnrmaResult <- cnrma(celFiles, cdfName="genomewidesnp6")

 		save(cnrmaResult, file=file.path(outdir, "cnrmaResult.rda"))

 	}

 }

@@ -94,8 +94,9 @@ We require 6 items for copy number estimation:

 <<snpAndCnSummaries>>=

 A <- res$A

 B <- res$B

-calls <- crlmmResult$calls

-conf <- crlmmResult$conf

+genotypes <- crlmm:::calls(crlmmResult)##$calls

+conf <- confs(crlmmResult)##$conf

+gender <- crlmmResult$gender

 SNR <- crlmmResult$SNR

 NP <- cnrmaResult$NP

 gc()

@@ -107,21 +108,22 @@ A histogram of the signal to noise ratio for the HapMap samples:

 hist(SNR, xlab="SNR", main="")

 @ 

-We suggest excluding samples with a signal to noise ratio less than 5.

-As batch effects can be very large in the quantile-normalized

-intensities, we suggest adjusting for date or chemistry plate.  

 <<dates>>=

 dts <- celDates(celFiles)

 table(format(dts, "%d %b %Y"))

+SNRmin <- 5

 @ 

-

-Ideally, one would have 70+ files in a given batch. Here we make a table

-of date versus ancestry:

+We suggest excluding samples with a signal to noise ratio less than

+\Sexpr{SNRmin}.  As batch effects can be very large in the

+quantile-normalized intensities, we suggest adjusting for date or

+chemistry plate.  Ideally, one would have 70+ files in a given

+batch. Here we make a table of date versus ancestry:

 <<specifyBatch>>=

-sns <- colnames(calls)

+require(Biobase)

+sns <- sampleNames(crlmmResult)

 sns[1]

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

 table(plate)

@@ -142,7 +144,7 @@ if(!exists("env")){

 		computeCopynumber(chrom=CHR,

 				  A=A,

 				  B=B,

-				  calls=calls,

+				  calls=genotypes,

 				  conf=conf,

 				  NP=NP,

 				  plate=plate,

@@ -150,7 +152,7 @@ if(!exists("env")){

 				  DF.PRIOR=50, 

 				  MIN.OBS=1,

 				  SNR=SNR,

-				  SNRmin=5)	

+				  SNRmin=SNRmin)	

 		save(env, file=file.path(outdir, paste("env_chr", CHR, ".rda", sep="")))

 	} else {

 		load(file.path(outdir, paste("env_chr", CHR, ".rda", sep="")))

@@ -165,8 +167,9 @@ that all of the assay data is bound to the meta-data.

 <<oligoSnpSet>>=

 require(oligoClasses)

-data(snpProbes, package="genomewidesnp6Crlmm")

-data(cnProbes, package="genomewidesnp6Crlmm")

+path <- system.file("extdata", package="genomewidesnp6Crlmm")

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

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

 position <- snpProbes[match(env[["snps"]], rownames(snpProbes)), "position"]

 position.np <- cnProbes[match(rownames(env[["NP"]]), rownames(cnProbes)), "position"]

 CA <- env[["CA"]]

@@ -220,7 +223,7 @@ if(!exists("env")){

 			computeCopynumber(chrom=CHR,

 					  A=A,

 					  B=B,

-					  calls=calls,

+					  calls=genotypes,

 					  conf=conf,

 					  NP=NP,

 					  plate=plate,

@@ -332,8 +335,8 @@ i <- I[1]

 ##plate effects are small

 log2(phiA[i, ])

 log2(phiB[i, ])

-plot(log2(A[i, ]), log2(B[i, ]), pch=as.character(calls[i, ]), col="grey60", cex=0.9, ylim=ylim, xlim=ylim, xlab="A", ylab="B")

-points(log2(A[i, J]), log2(B[i, J]), col="black", pch=as.character(calls[i, J]))

+plot(log2(A[i, ]), log2(B[i, ]), pch=as.character(genotypes[i, ]), col="grey60", cex=0.9, ylim=ylim, xlim=ylim, xlab="A", ylab="B")

+points(log2(A[i, J]), log2(B[i, J]), col="black", pch=as.character(genotypes[i, J]))

 for(CT in 2){

 	if(CT == 2) ellipse.col <- "black" ##else ellipse.col <- "purple"

 	for(CA in 0:CT){

@@ -419,9 +422,9 @@ for(j in seq(along=index)){

 	J1 <- grep(plate[plate1[j]], plate) ##sample indices for this plate

 	J2 <- grep(plate[plate2[j]], plate) ##sample indices for this plate

 	ylim <- c(6.5,13)

-	plot(log2(A[i, ]), log2(B[i, ]), pch=as.character(calls[i, ]), col="grey60", cex=0.9, ylim=ylim, xlim=ylim, xlab="A", ylab="B")

-	points(log2(A[i, J1]), log2(B[i, J1]), col="brown", pch=as.character(calls[i, J1]))

-	points(log2(A[i, J2]), log2(B[i, J2]), col="blue", pch=as.character(calls[i, J2]))

+	plot(log2(A[i, ]), log2(B[i, ]), pch=as.character(genotypes[i, ]), col="grey60", cex=0.9, ylim=ylim, xlim=ylim, xlab="A", ylab="B")

+	points(log2(A[i, J1]), log2(B[i, J1]), col="brown", pch=as.character(genotypes[i, J1]))

+	points(log2(A[i, J2]), log2(B[i, J2]), col="blue", pch=as.character(genotypes[i, J2]))

 	CT <- 2

 	P <- c(plate1[j], plate2[j])

 	for(p in seq(along=P)){

@@ -459,11 +462,6 @@ for(i in seq(along=index)){

 }

 @ 

-%\begin{figure}

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

-%  \caption{The prediction regions for copy number 2 shift when there

-%    is a large plate effect.}

-%\end{figure}

 \section{Session information}

 <<>>=

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

@@ -6,12 +6,12 @@

   to a HapMap reference distribution.

 }

 \usage{

-cnrma(filenames, sns, cdfName, seed = 1, verbose=FALSE)

+cnrma(filenames, cdfName="genomewidesnp6", sns, seed = 1, verbose=FALSE)

 }

 \arguments{

   \item{filenames}{filenames with complete path}

+  \item{cdfName}{Only 'genomewidesnp6' allowed}

   \item{sns}{the sample names. If missing, the basename of the filenames is used.}

-  \item{cdfName}{Should be 'genomewidesnp6Crlmm'}

   \item{seed}{seed for sampling the intensities to calculate skewness}

   \item{verbose}{logical}

 }

@@ -22,13 +22,11 @@ cnrma(filenames, sns, cdfName, seed = 1, verbose=FALSE)

 \author{Rob Scharpf}

 \note{Not tested}

 \examples{

-	\dontrun{

-		library(genomewidesnp6Crlmm)

-		library(hapmapsnp6)

-		path <- system.file("celFiles", package="hapmapsnp6")

-		celFiles <- list.celfiles(path, full.names=TRUE)

-		cnrmaResult <- cnrma(celFiles)

-	}

+	library(genomewidesnp6Crlmm)

+	library(hapmapsnp6)

+	path <- system.file("celFiles", package="hapmapsnp6")

+	celFiles <- list.celfiles(path, full.names=TRUE)

+	cnrmaResult <- cnrma(celFiles, cdfName="genomewidesnp6")

 }

 \keyword{robust}

@@ -13,7 +13,6 @@ computeCopynumber(chrom, A, B, calls, conf, NP, plate, MIN.OBS=1,

 envir, P, DF.PRIOR = 50, CONF.THR = 0.99, bias.adj=FALSE,

 priorProb, gender=NULL, SNR, SNRmin, seed=123, verbose=TRUE, ...)

 }

-%- maybe also 'usage' for other objects documented here.

 \arguments{

   \item{chrom}{Chromosome (an integer).  Use 23 for X and 24 for Y.}

   \item{A}{The A allele intensities from \code{snprma}}

@@ -52,22 +51,58 @@ priorProb, gender=NULL, SNR, SNRmin, seed=123, verbose=TRUE, ...)

 \details{

-  This function essentially transforms intensities from the raw scale to

-  a copy number scale.  We assume that the median within-SNP intensity

-  across samples is 2.  We make no assumption about the chromosomal copy

-  number. This function is mostly useful for detecting rare variants

-  (e.g., variants that would not affect the SNP-specific quantile-based

-  estimators of location and scale).  A correction for more common

-  variants is coming, as well as improved estimates of the uncertainty.

+  This function transforms the intensities to a copy number scale.  We

+  assume that the median within-SNP intensity across samples is 2.  We

+  make no assumption about the chromosomal copy number. This function is

+  useful for detecting rare variants (e.g., variants that would not

+  affect the SNP-specific quantile-based estimators of location and

+  scale).  A correction for more common variants is coming, as well as

+  improved estimates of the uncertainty.

 }

 \value{

-  For now, nothing is returned. All objects created by this function are

-  stored in the environment passed to this function.  See the copynumber

-  vignette for details on accessing these objects.

-  

+All objects created by this function are stored in the environment

+  passed to this function.  In addition, each of the elements are

+  specific to the chromosome(s) specified by the argument \code{chrom}.

+  For instance the element \code{A} is the matrix of quantile-normalized

+  intensities for the A-allele on chromosome(s) \code{chrom}.  The

+  element of this environment are as follows

+

+  \item{A}{Matrix of quantile-normalized intensities for the A-allele}

+  \item{B}{Matrix of quantile-normalized intensities for the A-allele}

+  \item{CA}{Copy number estimate for the A-allele (x 100)} 

+  \item{CB}{Copy number estimate for the B-allele (x 100)}

+  \item{calls}{CRLMM genotype calls (AA=1, AB=2, BB=3)}

+  \item{chrom}{Integer(s) indicating the chromosome(s)} 

+  \item{cnvs}{Names of the nonpolymorphic probes.  These are the

+  rownames of \code{NP} and \code{CT}.}

+  \item{conf}{CRLMM confidence scores for the genotypes: 'round(-1000*log2(1-p))'}

+  \item{corr}{Correlation of the A and B alleles for genotypes AB}

+  \item{corrA.BB}{Correlation of A and B alleles for genotypes BB}

+  \item{corrB.AA}{Correlation of A and B alleles for genotypes AA}

+  \item{CT}{Copy number estimates for nonpolymorphic probe.  See

+  \code{cnvs} for the rownames.}

+  \item{CT.sds}{Standard deviation estimates for \code{CT}}

+  \item{npflags}{Flags for the nonpolymorphic probes.}

+  \item{Ns}{The number of observations for each genotype/plate}

+  \item{nuA}{Background/cross-hyb for the A allele (plate- and locus-specific)}

+  \item{nuB}{Background/cross-hyb for the B allele (plate- and locus-specific)}

+  \item{nuT}{Background for the nonpolymorphic probes (plate- and locus-specific)}

+  \item{phiA}{Slope for the A allele (plate- and locus-specific)}

+  \item{phiB}{Slope for the B allele (plate- and locus-specific)}

+  \item{phiT}{Slope for the nonpolymorphic probes (plate- and locus-specific)}

+  \item{plate}{Factor indicating batch (same length as number of cel files)}

+  \item{sig2A}{Variance estimate for the A-allele signal (plate- and locus-specific)}

+  \item{sig2B}{Variance estimate for the B-allele signal (plate- and locus-specific)}

+  \item{sig2T}{Variance estimate for the nonpolymorphic signal (plate- and locus-specific)}

+  \item{snpflags}{Flags for polymorphic probes}

+  \item{snps}{Rownames for \code{A}, \code{B}, \code{CA}, \code{CB}, ...}

+  \item{sns}{ Sample names -- the column names for \code{A}, \code{B}, ...}

+  \item{steps}{Steps completed. For internal use.}

+  \item{tau2A}{Variance estimate for the B-allele background/cross-hyb (plate- and locus-specific)}

+  \item{tau2B}{Variance estimate for the B-allele background/cross-hyb (plate- and locus-specific)}

 }

 \references{Nothing yet.}

 \author{Rob Scharpf}