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

 .auto*

 test/*

+*.o

+*.so

 myCrlmmGT2.Rnw

 ffobjs*

 auto*

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

+Version: 1.9.20

 Date: 2010-12-10

 Author: Benilton S Carvalho <Benilton.Carvalho@cancer.org.uk>, Robert Scharpf <rscharpf@jhsph.edu>, Matt Ritchie <mritchie@wehi.edu.au>, Ingo Ruczinski <iruczins@jhsph.edu>, Rafael A Irizarry

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

@@ -70,4 +70,4 @@ export(crlmm,

 export(constructIlluminaCNSet)

 export(totalCopynumber, rawCopynumber)

 exportMethods(A, B, nuA, nuB, phiA, phiB, corr, tau2, Ns, medians, mads)

-

+export(constructInf, preprocessInf, genotypeInf)

@@ -100,21 +100,21 @@ construct <- function(filenames,

 }

 genotype <- function(filenames,

-		       cdfName,

-		       batch,

-		       mixtureSampleSize=10^5,

-		       eps=0.1,

-		       verbose=TRUE,

-		       seed=1,

-		       sns,

-		       probs=rep(1/3, 3),

-		       DF=6,

-		       SNRMin=5,

-		       recallMin=10,

-		       recallRegMin=1000,

+		     cdfName,

+		     batch,

+		     mixtureSampleSize=10^5,

+		     eps=0.1,

+		     verbose=TRUE,

+		     seed=1,

+		     sns,

+		     probs=rep(1/3, 3),

+		     DF=6,

+		     SNRMin=5,

+		     recallMin=10,

+		     recallRegMin=1000,

 		     gender=NULL,

-		       returnParams=TRUE,

-		       badSNP=0.7){

+		     returnParams=TRUE,

+		     badSNP=0.7){

 	is.lds <- ifelse(isPackageLoaded("ff"), TRUE, FALSE)

 	if(!is.lds) stop("this function now requires that the ff package be loaded")

 	if(missing(cdfName)) stop("must specify cdfName")

@@ -498,6 +498,7 @@ fit.lm1 <- function(strata,

 		    MIN.PHI,

 		    verbose, is.lds,

 		    CHR.X, ...){

+	open(object$gender)

 	if(is.lds) {physical <- get("physical"); open(object)}

 	if(verbose) message("      Probe stratum ", strata, " of ", length(index.list))

 	snps <- index.list[[strata]]

@@ -606,6 +607,7 @@ fit.lm2 <- function(strata,

 		    MIN.PHI,

 		    verbose, is.lds, CHR.X, ...){

 	if(is.lds) {physical <- get("physical"); open(object)}

+	open(object$gender)

 	if(verbose) message("      Probe stratum ", strata, " of ", length(index.list))

 	marker.index <- index.list[[strata]]

 	batches <- split(seq_along(batch(object)), as.character(batch(object)))

@@ -663,7 +665,9 @@ summarizeMaleXNps <- function(marker.index,

 	nr <- length(marker.index)

 	nc <- length(batchNames(object))

 	NN.Mlist <- imputed.medianA <- imputed.medianB <- shrink.madA <- shrink.madB <- vector("list", nc)

-	gender <- object$gender

+	open(object$gender)

+	gender <- object$gender[]

+	open(A(object))

 	AA <- as.matrix(A(object)[marker.index, gender==1])

 	madA.AA <- medianA.AA <- matrix(NA, nr, nc)

 	colnames(madA.AA) <- colnames(medianA.AA) <- batchNames(object)

@@ -697,10 +701,11 @@ summarizeXGenotypes <- function(marker.index,

 				DF.PRIOR,

 				gender="male", ...){

 	I <- unlist(batches)

+	open(object$gender)

 	if(gender == "male"){

-		gender.index <- which(object$gender == 1)

+		gender.index <- which(object$gender[] == 1)

 	} else {

-		gender.index <- which(object$gender == 2)

+		gender.index <- which(object$gender[] == 2)

 	}

 	batches <- lapply(batches, function(x, gender.index) intersect(x, gender.index), gender.index)

 	batch.names <- names(batches)

@@ -709,6 +714,10 @@ summarizeXGenotypes <- function(marker.index,

 	nc <- length(batch.index)

 	NN.Mlist <- medianA <- medianB <- shrink.madA <- shrink.madB <- vector("list", nc)

 	names(NN.Mlist) <- names(medianA) <- names(medianB) <- names(shrink.madA) <- names(shrink.madB) <- batch.names

+	open(calls(object))

+	open(snpCallProbability(object))

+	open(A(object))

+	open(B(object))

 	GG <- as.matrix(calls(object)[marker.index, ])

 	CP <- as.matrix(snpCallProbability(object)[marker.index, ])

 	AA <- as.matrix(A(object)[marker.index, ])

@@ -819,6 +828,10 @@ summarizeXGenotypes <- function(marker.index,

 		medianA[[k]] <- res[[1]]

 		medianB[[k]] <- res[[2]]

 	}

+	close(calls(object))

+	close(snpCallProbability(object))

+	close(A(object))

+	close(B(object))

 	return(list(madA=shrink.madA,

 		    madB=shrink.madB,

 		    NN=NN.Mlist,

@@ -841,7 +854,8 @@ fit.lm3 <- function(strata,

 		    verbose, is.lds, CHR.X, ...){

 	##if(is.lds) {physical <- get("physical"); open(object)}

 	if(verbose) message("      Probe stratum ", strata, " of ", length(index.list))

-	gender <- object$gender

+	open(object$gender)

+	gender <- object$gender[]

 	enough.males <- sum(gender==1) >= MIN.SAMPLES

 	enough.females <- sum(gender==2) >= MIN.SAMPLES

 	if(!enough.males & !enough.females){

@@ -1011,7 +1025,8 @@ fit.lm4 <- function(strata,

 		    verbose, is.lds=TRUE, ...){

 	##if(is.lds) {physical <- get("physical"); open(object)}

 	## exclude batches that have fewer than MIN.SAMPLES

-	gender <- object$gender

+	open(object$gender)

+	gender <- object$gender[]

 	enough.males <- sum(gender==1) >= MIN.SAMPLES

 	enough.females <- sum(gender==2) >= MIN.SAMPLES

 	if(!enough.males & !enough.females){

@@ -1434,7 +1449,7 @@ shrinkSummary <- function(object,

 			  DF.PRIOR=50,

 			  verbose=TRUE,

 			  marker.index,

-			  is.lds){

+			  is.lds=TRUE){

 	stopifnot(type[[1]] == "SNP")

 	CHR.X <- FALSE ## this is no longer needed

 	if(missing(marker.index)){

@@ -1469,7 +1484,7 @@ shrinkSummary <- function(object,

 			      DF.PRIOR=DF.PRIOR,

 			      is.lds=is.lds)

 	}

-	return(object)

+	TRUE

 }

 genotypeSummary <- function(object,

@@ -1522,7 +1537,7 @@ genotypeSummary <- function(object,

 			      is.lds=is.lds,

 			      CHR.X=CHR.X)

 	}

-	return(object)

+	TRUE

 }

 whichMarkers <- function(type, is.snp, is.autosome, is.annotated, is.X){

@@ -1552,6 +1567,8 @@ summarizeNps <- function(strata, index.list, object, batchSize,

 	nc <- length(batchnames)

 	N.AA <- medianA.AA <- madA.AA <- tau2A.AA <- matrix(NA, nr, nc)

 	is.illumina <- whichPlatform(annotation(object)) == "illumina"

+	open(A(object))

+	open(object$gender)

 	AA <- as.matrix(A(object)[index, ])

 	if(is.illumina){

 		BB <- as.matrix(B(object)[index, ])

@@ -1593,7 +1610,6 @@ summarizeSnps <- function(strata,

 ##		physical <- get("physical")

 ##		open(object)

 ##	}

-	open(object)

 	if(verbose) message("      Probe stratum ", strata, " of ", length(index.list))

 	index <- index.list[[strata]]

 	batches <- split(seq_along(batch(object)), as.character(batch(object)))

@@ -1789,6 +1805,7 @@ crlmmCopynumber <- function(object,

 	is.lds <- is(calls(object), "ffdf") | is(calls(object), "ff_matrix")

 	if(is.lds) stopifnot(isPackageLoaded("ff"))

 	samplesPerBatch <- table(as.character(batch(object)))

+	open(object)

 	if(any(samplesPerBatch < MIN.SAMPLES)){

 		tmp <- paste(names(samplesPerBatch)[samplesPerBatch < MIN.SAMPLES], collapse=", ")

 		message("The following batches have fewer than ", MIN.SAMPLES, " samples: ",  tmp)

@@ -1810,25 +1827,25 @@ crlmmCopynumber <- function(object,

 		##if(verbose) message(paste("Computing summary statistics for ", mylabel(marker.type), " genotype clusters for each batch")

 		marker.index <- whichMarkers(marker.type, is.snp,

 					     is.autosome, is.annotated, is.X)

-		object <- genotypeSummary(object=object,

-					  GT.CONF.THR=GT.CONF.THR,

-					  type=marker.type,

-					  verbose=verbose,

-					  marker.index=marker.index,

-					  is.lds=is.lds)

+		genotypeSummary(object=object,

+				GT.CONF.THR=GT.CONF.THR,

+				type=marker.type,

+				verbose=verbose,

+				marker.index=marker.index,

+				is.lds=is.lds)

 	}

 	if(verbose) message("Imputing unobserved genotype medians and shrinking the variances (within-batch, across loci) ")##SNPs only

 	if("SNP" %in% type){

 		marker.index <- whichMarkers("SNP", is.snp,

 					     is.autosome, is.annotated, is.X)

-		object <- shrinkSummary(object=object,

-					MIN.OBS=MIN.OBS,

-					MIN.SAMPLES=MIN.SAMPLES,

-					DF.PRIOR=DF.PRIOR,

-					type="SNP",

-					verbose=verbose,

-					marker.index=marker.index,

-					is.lds=is.lds)

+		shrinkSummary(object=object,

+			      MIN.OBS=MIN.OBS,

+			      MIN.SAMPLES=MIN.SAMPLES,

+			      DF.PRIOR=DF.PRIOR,

+			      type="SNP",

+			      verbose=verbose,

+			      marker.index=marker.index,

+			      is.lds=is.lds)

 	}

 	if(verbose) message("Estimating parameters for copy number")##SNPs only

 	for(i in seq_along(type)){

@@ -1852,6 +1869,7 @@ crlmmCopynumber <- function(object,

 					       is.lds=is.lds,

 					       CHR.X=CHR.X)

 	}

+	close(object)

 	TRUE

 }

 crlmmCopynumber2 <- function(){

@@ -1008,16 +1008,14 @@ getProtocolData.Illumina = function(filenames, sep="_", fileExt="Grn.idat", verb

 }

-construct.Illumina <- function(sampleSheet=NULL,

+constructInf <- function(sampleSheet=NULL,

 			       arrayNames=NULL,

-			       ids=NULL,

 			       path=".",

 			       arrayInfoColNames=list(barcode="SentrixBarcode_A", position="SentrixPosition_A"),

 			       highDensity=FALSE,

 			       sep="_",

 			       fileExt=list(green="Grn.idat", red="Red.idat"),

 			       cdfName,

-			       copynumber=TRUE,

 			       verbose=FALSE,

 			       batch, #fns,

 			       saveDate=TRUE) { #, outdir="."){

@@ -1082,14 +1080,8 @@ construct.Illumina <- function(sampleSheet=NULL,

 	if(!all(file.exists(redidats))) stop("Missing some of the *Red.idat files")

 	if(verbose) message("Initializing container for genotyping and copy number estimation")

-	featureData = getFeatureData(cdfName, copynumber=copynumber)

-#	if(!missing(fns)){

-#		index = match(fns, featureNames(featureData))

-#		if(all(is.na(index))) stop("fns not in featureNames")

-#		featureData = featureData[index, ]

-#	}

+	featureData = getFeatureData(cdfName, copynumber=TRUE)

 	nr = nrow(featureData); nc = narrays

-#        ldPath(outdir)

 	cnSet <- new("CNSet",

 		     alleleA=initializeBigMatrix(name="A", nr, nc),

 		     alleleB=initializeBigMatrix(name="B", nr, nc),

@@ -1114,13 +1106,12 @@ construct.Illumina <- function(sampleSheet=NULL,

 	cnSet$SNR = initializeBigVector("crlmmSNR-", ncol(cnSet), "double")

 	cnSet$SKW = initializeBigVector("crlmmSKW-", ncol(cnSet), "double")

 	sampleNames(cnSet) <- basename(sampleNames(cnSet))

-	##pd = data.frame(matrix(NA, nc, 3), row.names=sampleNames(cnSet))

-	##colnames(pd)=c("SKW", "SNR", "gender")

-	##phenoData(cnSet) = new("AnnotatedDataFrame", data=pd)

 	return(cnSet)

 }

+construct.Illumina <- constructInf

-preprocess <- function(cnSet,

+preprocessInf <- function(cnSet,

+			  sampleSheet=NULL,

 		       arrayNames=NULL,

 		       ids=NULL,

 		       path=".",

@@ -1176,16 +1167,17 @@ preprocess <- function(cnSet,

 		 neededPkgs=c("crlmm", pkgname)) # outdir=outdir,

 	return(mixtureParams)

 }

+preprocess <- preprocessInf

-genotypeRS <- function(cnSet, mixtureParams, probs=rep(1/3,3),

-		       SNRMin=5,

-		       recallMin=10,

-		       recallRegMin=1000,

-		       verbose=TRUE,

-		       returnParams=TRUE,

-		       badSNP=0.7,

-		       gender=NULL,

-		       DF=6){

+genotypeInf <- function(cnSet, mixtureParams, probs=rep(1/3,3),

+			SNRMin=5,

+			recallMin=10,

+			recallRegMin=1000,

+			verbose=TRUE,

+			returnParams=TRUE,

+			badSNP=0.7,

+			gender=NULL,

+			DF=6){

 	is.snp = isSnp(cnSet)

 	snp.index = which(is.snp)

 	narrays = ncol(cnSet)

@@ -1256,91 +1248,49 @@ genotype.Illumina <- function(sampleSheet=NULL,

 	stopifnot(is.lds)

 	if(missing(cdfName)) stop("must specify cdfName")

 	if(!isValidCdfName(cdfName)) stop("cdfName not valid.  see validCdfNames")

-        pkgname = getCrlmmAnnotationName(cdfName)

-	callSet <- construct.Illumina(sampleSheet=sampleSheet, arrayNames=arrayNames,

-				      ids=ids, path=path, arrayInfoColNames=arrayInfoColNames,

-				      highDensity=highDensity, sep=sep, fileExt=fileExt,

-				      cdfName=cdfName, copynumber=copynumber, verbose=verbose, batch=batch, # fns=fns,

-				      saveDate=saveDate) #, outdir=outdir)

-	## Basic checks

-	##check that mixtureParams provided

-	if(missing(mixtureParams))

-		stop("preprocess is FALSE but mixtureParams were not provided")

-	if(ncol(callSet) != length(arrayNames))

-		stop("callSet provided but number of samples is not the same as the length of arrayNames")

-	snr <- callSet$SNR[]

-	if(any(is.na(snr))) stop("missing values in callSet$SNR")

-	if(missing(sns)) sns = basename(sampleNames(callSet))

-	open(A(callSet))

-	open(B(callSet))

- 	is.snp = isSnp(callSet)

-	snp.index = which(is.snp)

-	narrays = ncol(callSet)

-	sampleBatches <- splitIndicesByLength(seq(length=ncol(callSet)), ocSamples())

-	mixtureParams = initializeBigMatrix("crlmmMixt-", 4, narrays, "double")

-	SNR = initializeBigVector("crlmmSNR-", narrays, "double")

-	SKW = initializeBigVector("crlmmSKW-", narrays, "double")

-	ocLapply(seq_along(sampleBatches), processIDAT, sampleBatches=sampleBatches,

-		 sampleSheet=sampleSheet, arrayNames=arrayNames,

-		 ids=ids, path=path, arrayInfoColNames=arrayInfoColNames, highDensity=highDensity,

-		 sep=sep, fileExt=fileExt, saveDate=saveDate, verbose=verbose, mixtureSampleSize=mixtureSampleSize,

-		 fitMixture=fitMixture, eps=eps, seed=seed, cdfName=cdfName, sns=sns, stripNorm=stripNorm,

-		 useTarget=useTarget, A=A(callSet), B=B(callSet), SKW=SKW, SNR=SNR,

-		 mixtureParams=mixtureParams, is.snp=is.snp, neededPkgs=c("crlmm", pkgname)) # outdir=outdir,

-	open(SKW)

-	open(SNR)

-	pData(callSet)$SKW = SKW

-	pData(callSet)$SNR = SNR

-	save(callSet, file=file.path(ldPath(), "callSet.rda"))

-	close(SNR)

-	close(SKW)

-	open(A(callSet))

-	open(B(callSet))

-	tmpA = initializeBigMatrix(name="tmpA", length(snp.index), narrays)

-	tmpB = initializeBigMatrix(name="tmpB", length(snp.index), narrays)

-	for(j in seq(length=ncol(callSet))){

-		tmpA[, j] <- as.integer(A(callSet)[snp.index, j])

-		tmpB[, j] <- as.integer(B(callSet)[snp.index, j])

-	}

-	close(A(callSet))

-	close(B(callSet))

-	close(tmpA)

-	close(tmpB)

-	save(tmpA, file=file.path(ldPath(), "tmpA.rda"))

-	save(tmpB, file=file.path(ldPath(), "tmpB.rda"))

-	save(SNR, file=file.path(ldPath(), "SNR.rda"))

-	save(SKW, file=file.path(ldPath(), "SKW.rda"))

-	save(mixtureParams, file=file.path(ldPath(), "mixtureParams.rda"))

-	message("Begin genotyping")

-	tmp <- crlmmGT2(A=tmpA,

-			B=tmpB,

-			SNR=SNR,

-			mixtureParams=mixtureParams,

-			cdfName=cdfName,

-			row.names=NULL,

-			col.names=sampleNames(callSet),

-			probs=probs,

-			DF=DF,

-			SNRMin=SNRMin,

-			recallMin=recallMin,

-			recallRegMin=recallRegMin,

-			gender=gender,

-			verbose=verbose,

-			returnParams=returnParams,

-			badSNP=badSNP)

-	save(tmp, file=file.path(ldPath(), "tmp.rda"))

-	if(verbose) message("Genotyping finished.  Updating container with genotype calls and confidence scores.")

-	open(tmpA); open(tmpB)

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

-		snpCall(callSet)[snp.index, j] <- as.integer(tmpA[, j])

-		snpCallProbability(callSet)[snp.index, j] <- as.integer(tmpB[, j])

-	}

-	close(tmpA); close(tmpB)

-	delete(tmpA); delete(tmpB);

-	callSet$gender = tmp$gender

-	rm(tmp)

-	close(callSet)

-	return(callSet)

+	stopifnot(!missing(batch))

+        pkgname <- getCrlmmAnnotationName(cdfName)

+	message("Instantiate CNSet container.")

+	cnSet <- constructInf(sampleSheet=sampleSheet,

+				    arrayNames=arrayNames,

+				    path=path,

+				    arrayInfoColNames=arrayInfoColNames,

+				    highDensity=highDensity,

+				    sep=sep,

+				    fileExt=fileExt,

+				    cdfName=cdfName,

+				    verbose=verbose,

+				    batch=batch,

+				    saveDate=saveDate)

+	mixtureParams <- preprocessInf(cnSet=cnSet,

+				    sampleSheet=sampleSheet,

+				    arrayNames=arrayNames,

+				    ids=ids,

+				    path=path,

+				    arrayInfoColNames=arrayInfoColNames,

+				    highDensity=highDensity,

+				    sep=sep,

+				    fileExt=fileExt,

+				    saveDate=saveDate,

+				    stripNorm=stripNorm,

+				    useTarget=useTarget,

+				    mixtureSampleSize=mixtureSampleSize,

+				    fitMixture=fitMixture,

+				    eps=eps,

+				    verbose=verbose,

+				    seed=seed)

+	message("Preprocessing complete.  Begin genotyping...")

+	genotypeInf(cnSet=cnSet, mixtureParams=mixtureParams,

+		   probs=probs,

+		   SNRMin=SNRMin,

+		   recallMin=recallMin,

+		   recallRegMin=recallRegMin,

+		   verbose=verbose,

+		   returnParams=returnParams,

+		   badSNP=badSNP,

+		   gender=gender,

+		   DF=DF)

+	return(cnSet)

 }

@@ -55,3 +55,44 @@ linesCNSet <- function(x, y, batch, copynumber, x.axis="A", grid=FALSE, ...){

 		}

 	}

 }

+

+

+cnPanel <- function(x, y, ..., pch.cols, gt, cbs.segs, hmm.segs=NULL, shades, subscripts, add.ideogram=TRUE){

+	##if(panel.number() == 2) browser()

+	add.ideogram <- add.ideogram[[panel.number()]]

+	##cbs.segs <- cbs.segs[[panel.number()]]

+	draw.hmm.states <- ifelse(panel.number() <= length(hmm.segs), TRUE, FALSE)

+	panel.grid(h=6, v=10)

+	which.hom <- which(gt[subscripts] == 1 | gt[subscripts]==3)

+	which.het <- which(gt[subscripts] == 2)

+	panel.xyplot(x, y, col="grey60", ...)

+	lpoints(x[which.hom], y[which.hom], col=pch.cols[1], ...)

+	lpoints(x[which.het], y[which.het], col=pch.cols[2], ...)

+	lsegments(x0=start(cbs.segs)/1e6, x1=end(cbs.segs)/1e6,

+		  y0=cbs.segs$seg.mean,

+		  y1=cbs.segs$seg.mean, ...)

+	if(draw.hmm.states){

+		hmm.segs <- hmm.segs[order(width(hmm.segs), decreasing=TRUE), ]

+		lrect(xleft=start(hmm.segs)/1e6,

+		      xright=end(hmm.segs)/1e6,

+		      ybottom=-0.4,

+		      ytop=0,

+		      border=shades[hmm.segs$state],

+		      col=shades[hmm.segs$state])

+	}

+	ltext(130, 5, paste("MAD:", round(mad(y, na.rm=TRUE), 2)))

+	if(add.ideogram){

+		pathto <- system.file("hg18", package="SNPchip")

+		cytoband <- read.table(file.path(pathto, "cytoBand.txt"), as.is=TRUE)

+		cytoband$V2 <- cytoband$V2/1e6

+		cytoband$V3 <- cytoband$V3/1e6

+		colnames(cytoband) <- c("chrom", "start", "end", "name", "gieStain")

+		plotCytoband(unique(hmm.segs$chrom),

+			     cytoband=cytoband,

+			     cytoband.ycoords=c(5.6, 5.9),

+			     new=FALSE,

+			     label.cytoband=FALSE,

+			     build="hg18",

+			     use.lattice=TRUE)

+	}

+}

@@ -15,6 +15,7 @@

 \newcommand{\Rclass}[1]{{\textit{#1}}}

 \newcommand{\oligo}{\Rpackage{oligo }}

 \newcommand{\R}{\textsf{R}}

+\newcommand{\ff}{\Rpackage{ff}}

 \begin{document}

 \title{Using \crlmm{} for copy number estimation and genotype calling

@@ -27,619 +28,291 @@

 \begin{abstract}

-  This vignette illustrates the steps necessary for obtaining

-  marker-level estimates of allele-specific copy number from the raw

-  Illumina IDAT files.  Hidden markov models or segmentation

-  algorithms can be applied to the marker-level estimates.  Examples

-  of both are illustrated.

+  This vignette illustrates the steps required prior to copy number

+  analysis for Infinium platforms.  Specifically, we require

+  construction of a container to store processed forms of the raw

+  data, preprocessing to normalize the arrays, and genotyping using

+  the CRLMM algorithm.  After completing these steps, users can refer

+  to the copy number vignette.

 \end{abstract}

-\section{About this vignette}

+\section{Infrastructure}

+

 <<crlmm, results=hide, echo=FALSE>>=

 library(crlmm)

+options(width=70)

+options(continue=" ")

 @

-%The vignette for copy number estimation for illumina platforms is

-%under development.  Please use the latest stable release of crlmm.

-

-%\end{document}

+\textbf{Supported platforms:} The supported Infinium platforms are

+those for which a corresponding annotation package is available.  The

+annotation packages contain information on the markers, such as

+physical position and chromosome, as well as pre-computed parameters

+estimated from HapMap used during the preprocessing and genotyping

+steps. Currently supported Infinium platforms are listed in the

+following code chunk.

-Allele-specific copy number estimation in the crlmm package is

-available for several Illumina platforms.  As described in the

-\texttt{copynumber} vignette, copy number estimation in \crlmm{} works

-best when there are a sufficient number of samples such that AA, AB,

-and BB genotypes are observed at most loci. For small studies (e.g.,

-fewer than 50 samples), there will be a large number of SNPs that are

-monomorphic. For monomorphic SNPs, the estimation problem becomes more

-difficult and alternative strategies that estimate the relative total

-copy number may be preferable.  In addition to installing \crlmm{},

-one must also install the appropriate annotation package for the

-Illumina platform.  In the following code, we list the platforms for

-which annotation packages are currently available. Next we create a

-directory where output files will be stored and indicate the directory

-that contains the IDAT files that will be used in our analysis.

-

-

-<<setup, echo=FALSE, results=hide>>=

-options(width=70)

-options(continue=" ")

-library(cacheSweave)

+<<supportedPlatforms>>=

+pkgs <- annotationPackages()

+crlmm.pkgs <- pkgs[grep("Crlmm", pkgs)]

+crlmm.pkgs[grep("human", crlmm.pkgs)]

 @

-<<libraries>>=

+\textbf{Large data:} In order to reduce \crlmm{}'s memory footprint,

+we require the \Rpackage{ff} for copy number estimation.  The

+\Rpackage{ff} package provides infrastructure for accessing and

+writing data to disk instead of keeping data in memory.  Each element

+of the \verb+assayData+ and \verb+batchStatistics+ slot of the

+\Rclass{CNSet} class are \Robject{ff} objects.  \Robject{ff} objects

+in the \R{} workspace contain pointers to several files with the

+\verb+.ff+ extension.  The location of where the data is stored on

+disk is specified by use of the \Rfunction{ldPath} function.

+

+<<ldpath,results=hide>>=

 library(ff)

-pkgs <- annotationPackages()

-pkgs[grep("Crlmm", pkgs)]

 if(getRversion() < "2.13.0"){

 	rpath <- getRversion()

 } else rpath <- "trunk"

 outdir <- paste("/thumper/ctsa/snpmicroarray/rs/ProcessedData/crlmm/", rpath, "/illumina_vignette", sep="")

+ldPath(outdir)

 dir.create(outdir, recursive=TRUE, showWarnings=FALSE)

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

-setCacheDir(outdir)

 @

-%Options for controlling RAM with the \Rpackage{ff} package.

-

-<<ldOptions>>=

-ldPath(outdir)

+Users should not move or rename this directory.  If only output files

+are stored in \verb+outdir+, one can either remove the entire

+directory prior to rerunning the analysis or all of the '.ff' files.

+Otherwise, one would accumulate a large number of '.ff' files on disk

+that are no longer in use.

+

+As none of the functions for preprocessing, genotyping, or copy number

+estimation simultaneously require all samples and all probes in

+memory, memory-usage by \crlmm{} can be fine tuned by reading in and

+processing subsets of the markers and/or samples. The functions

+\Rfunction{ocSamples} and \Rfunction{ocProbesets} in the

+\Rpackage{oligoClasses} package can be used to declare how many

+markers and samples to read at once.  In general, specifying smaller

+values should reduce the RAM required for a particular job, but would

+be expected to have an increased run-time. In the following

+code-chunk, we declare that \crlmm{} should process 150,000 markers at

+a time (when possible) and 500 samples at a time.  (As our dataset in

+this vignette only contains 43 samples, the \Rfunction{ocSamples}

+option would not have any effect.)  One can view the current settings

+for these commands, by typing the functions without an argument.

+

+<<ram>>=

+ocSamples()

+ocProbesets()

 ocProbesets(150e3)

-ocSamples(200)

+ocSamples(500)

+ocSamples()

+ocProbesets()

 @

-This vignette was created using Illumina IDAT files that are located

-in a specific directory on my computer (\Robject{pathToCels}).  Long

-computations are saved in the output directory \Robject{outdir}.

-Towards this end, we make repeated use of the \Rfunction{checkExists}

-simply checks that a variable is in the workspace.  If not, the

-variable is loaded from the indicated path.  If the file (with '.rda'

-extension) does not exist in this path, the function indicated by the

-.FUN argument is executed.  Alternatively, if \Robject{.load.it} is

-\Robject{FALSE} the function will be executed regardless of whether

-the file exists.  Users should see the \Rpackage{weaver} package for a

-more 'correct' approach for cacheing long computations. The following

-code chunk checks that that the variables specific to the directory

-structure on my computer exist.  Users should modify the

-\Robject{outdir} and \Robject{datadir} variables as appropriate.

-

-<<checkSetup>>=

-if(!file.exists(outdir)) stop("Please specify valid directory for storing output")

-if(!file.exists(datadir)) stop("Please specify the correct path to the CEL files")

+<<cacheSweave, echo=FALSE, results=hide>>=

+library(cacheSweave)

+setCacheDir(outdir)

 @

-\section{Preprocessing Illumina IDAT files, genotyping, and copy number

-  estimation}

-To perform copy number analysis on the Illumina platform, several

-steps are required.  The first step is to read in the IDAT files and

-create a container for storing the red and green intensities. These

-intensities are quantile normalized in the function

-\Rfunction{crlmmIllumina}, and then genotyped using the crlmm

-algorithm.  Details on the crlmm genotyping algorithm are described

-elsewhere.  We will make use of the normalized intensities when we

-estimate copy number.  The object returned by the

-\Rfunction{genotype.Illumina} function is an instance of the

-\Robject{CNSet} class, a container for storing the genotype calls,

-genotype confidence scores, the normalized intensities for the A and B

-channels, and summary statistics on the batches.  The batch variable

-can be the date on which the array was scanned or the 96 well

-microarray chemistry plate on which the samples were processed. Both

-date and chemistry plate are useful surrogates for experimental

-factors that effect subgroups of the probe intensities.  (Quantile

-normalization does not remove batch effects.)  The genotype confidence

-scores are saved as an integer, and can be converted back to a $[0,

-1]$ probability scale by the transformation $round(-1000*log2(1-p))$.

+\textbf{Limitations:} There is no minimum number of samples required

+for preprocessing and genotyping.  However, for copy number estimation

+the \crlmm{} package currently requires at least 10 samples per batch.

+The parameter estimates for copy number and the corresponding

+estimates of raw copy number will tend to be more noisy for batches

+with small sample sizes (e.g., $<$ 50).  Chemistry plate or scan date

+are often useful surrogates for batch.

+

+\section{Initializing a container for storing processed data}

+

+This section will initialize a container for storing processed forms

+of the data, including the normalized intensities for the A and B

+alleles and the CRLMM genotype calls and confidence scores.  In

+addition, the container will store information on the markers

+(physical position, chromosome, and a SNP indicator), the batch, and

+the samples (e.g., gender).  To construct this container for Infinium

+platforms, several steps are required.

+

+We begin by specifying the path containing the raw IDAT files for a

+set of samples from the Infinium 370k platform.

+

+<libraries>>=

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

+@

+

+For Infinium platforms, an Illumina sample sheet containing

+information for reading the raw IDAT files is required. Please refer

+to the BeadStudio Genotyping guide, Appendix A, for additional

+information.  The following code reads in the samplesheet for the IDAT

+files on our local server. For our dataset, the file extensions are

+`Grn.dat' and `Red.idat'.  We store the complete path to the filename

+without the file extension in the \Robject{arrayNames} and check that

+all of the green and red IDAT files exists.

 <<samplesheet>>=

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

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

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

-grnfiles = all(file.exists(paste(arrayNames, "_Grn.idat", sep="")))

-redfiles = all(file.exists(paste(arrayNames, "_Red.idat", sep="")))

-load.it <- TRUE

-if(!load.it){

-	rmFiles <- list.files(outdir, pattern=".ff", full.names=TRUE)

-	unlink(rmFiles)

-}

-plate <- samplesheet$Plate

-table(plate)

-## too few samples to run the plates as separate batches

-batch <- rep("1", nrow(samplesheet))

+all(file.exists(paste(arrayNames, "_Grn.idat", sep="")))

+all(file.exists(paste(arrayNames, "_Red.idat", sep="")))

+arrayInfo <- list(barcode=NULL, position="SentrixPosition")

 @

-<<genotype,cache=TRUE>>=

-container <- genotype.Illumina(sampleSheet=samplesheet,

-			       path=dirname(arrayNames[1]),

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

-			       cdfName="human370v1c",

-			       batch=batch)

-@

+As discussed previously, all supported platforms have a corresponding

+annotation package.  The appropriate annotation package is specified

+by the platform identifier without the \verb+Crlmm+ postfix.

-<<copynumber,cache=TRUE>>=

-GT.CONF.THR <- 0.8

-cnSet <- crlmmCopynumber(container, GT.CONF.THR=GT.CONF.THR)

+<<cdfname>>=

+cdfName <- "human370v1c"

 @

-The \Robject{cnSet} returned by the \Rfunction{crlmmCopynumber}

-function contains all of the parameters used to compute

-allele-specific copy number.  In order to reduce I/O and speed

-computation, the allele-specific copy number estimates are not written

-to file nor saved in the \Robject{CNSet} object.  Rather,

-allele-specific estimates are computed on the fly by the functions

-\Rfunction{CA}, \Rfunction{CB}, or \Rfunction{totalCopynumber}.  The

-following codechunks provide a few examples.

-

-Copy number for polymorphic markers on chromosomes 1 - 22.

-

-<<ca>>=

-snp.index <- which(isSnp(cnSet) & !is.na(chromosome(cnSet)) & chromosome(cnSet) < 22)

-ca <- CA(cnSet, i=snp.index, j=1:5)

-cb <- CB(cnSet, i=snp.index, j=1:5)

-ct <- ca+cb

-@

+Next, we construct a character vector that specifies the batch for

+each of the \Sexpr{length(arrayNames)} arrays.  Here, we have a small

+dataset and process the samples in a single batch. Processing the

+samples as a single batch is generally reasonable if the samples were

+processed at similar times (e.g., within a few weeks).

-Alternatively, total copy number can be obtained by

-<<totalCopynumber.snps>>=

-ct2 <- totalCopynumber(cnSet, i=snp.index, j=1:5)

-stopifnot(all.equal(ct, ct2))

+<<batch>>=

+batch <- rep("1", nrow(samplesheet))

 @

-Missing values can arise at polymorphic when the confidence score of

-the genotype calls are below the threshold indicated by the threshold

-\texttt{GT.CONF.THR} in \Robject{crlmmCopynumber}. See

-\Robject{?crlmmCopynumber} for additional details.  In the following

-codechunk, we compute the number of samples that had confidence scores

-below \Sexpr{GT.CONF.THR} at loci for which ${\hat CA}$ and ${\hat CB}$

-are missing.

-

-<<NAs.snps>>=

-missing.copynumber <- which(rowSums(is.na(ct)) > 0)

-invisible(open(snpCallProbability(cnSet)))

-gt.confidence <- i2p(snpCallProbability(cnSet)[snp.index, ])

-n.below.threshold <- rowSums(gt.confidence < 0.95)

-unique(n.below.threshold[missing.copynumber])

-@

-\noindent For loci with missing copy number, the confidence scores

-were all below the threshold.

-

-At nonpolymorphic loci, either the \Rfunction{CA} or

-\Rfunction{totalCopynumber} functions can be used to obtain estimates

-of total copy number.

-

-<<nonpolymorphicAutosomes>>=

-marker.index <- which(!isSnp(cnSet) & chromosome(cnSet) < 23)

-ct <- CA(cnSet, i=marker.index, j=1:5)

-## all zeros

-stopifnot(all(CB(cnSet, i=marker.index, j=1:5) == 0))

-ct2 <- totalCopynumber(cnSet, i=marker.index, j=1:5)

-stopifnot(all.equal(ct, ct2))

+Finally, we initialize an object of class \Robject{CNSet} using the

+function \Rfunction{constructInf}.

+

+<<container,cache=TRUE>>=

+cnSet <- constructInf(sampleSheet=samplesheet,

+		      arrayNames=arrayNames,

+		      batch=batch,

+		      arrayInfoColNames=arrayInfo,

+		      cdfName=cdfName,

+		      verbose=TRUE,

+		      saveDate=TRUE)

 @

+A concise summary of the object's contents can be viewed with the

+\Rfunction{print} function.

-\begin{figure}[t!]

-  Nonpolymorphic markers on chromosome X:

-  \centering

-<<nonpolymorphicX, fig=TRUE, width=8, height=4>>=

-npx.index <- which(chromosome(cnSet)==23 & !isSnp(cnSet))

-M <- sample(which(cnSet$gender==1), 5)

-F <- sample(which(cnSet$gender==2), 5)

-cn.M <- CA(cnSet, i=npx.index, j=M)

-cn.F <- CA(cnSet, i=npx.index, j=F)

-boxplot(data.frame(cbind(cn.M, cn.F)), pch=".", col="grey60", outline=FALSE)

-@

-\caption{Copy number estimates for nonpolymorphic loci on chromosome

-  X (5 men, 5 women).  crlmm assumes that the median copy number

-  across samples at a given marker on X is 1 for men and 2 for women.}

-\end{figure}

-

-

-

-

-\begin{figure}[t!]

- \centering

-<<polymorphicX, fig=TRUE, width=8, height=4>>=

-## copy number estimates on X for SNPs are biased towards small values.

-X.markers <- which(isSnp(cnSet) & chromosome(cnSet) == 23)

-ca.M <- CA(cnSet, i=X.markers, j=M)

-cb.M <- CB(cnSet, i=X.markers, j=M)

-ca.F <- CA(cnSet, i=X.markers, j=F)

-cb.F <- CB(cnSet, i=X.markers, j=F)

-cn.M <- ca.M+cb.M

-cn.F <- ca.F+cb.F

-boxplot(data.frame(cbind(cn.M, cn.F)), pch=".", outline=FALSE, col="grey60")

-cn2 <- totalCopynumber(cnSet, i=X.markers, j=c(M, F))

-stopifnot(all.equal(cbind(cn.M, cn.F), cn2))

-@

-\caption{Copy number estimates for polymorphic markers on chromosome

-  X. }

-\end{figure}

-

-Often it is of interest to smooth the total copy number estimates (the

-sum of the allele-specific copy number) as a function of the physical

-position.  The following code chunk can be used to create an instance

-of \Robject{CopyNumberSet} that contains the CA + CB at polymorphic

-markers and 'CA' at nonpolymorphic markers.

-

-The \Robject{cnSet} returned by the \Rfunction{crlmmCopynumber}

-function contains all of the parameters used to compute

-allele-specific copy number.  In order to reduce I/O and speed

-computation, the allele-specific copy number estimates are not written

-to file nor saved in the \Robject{CNSet} object.  Rather, the

-estimates are computed on the fly when the user calls one of the

-helper functions that uses this information. For instance, often it is

-of interest to smooth the total copy number estimates (the sum of the

-allele-specific copy number) as a function of the physical position.

-The following code chunk can be used to create an instance of

-\Robject{CopyNumberSet} that contains the CA + CB at polymorphic

-markers and 'CA' at nonpolymorphic markers.

-