@@ -29,14 +29,13 @@

   This vignette provides an overview of the \Rclass{CNSet} class and a

   brief discussion of the underlying infrastructure for large data

-  support with the \Rpackage{ff} package.  This package instantiates

+  support with the \Rpackage{ff} package.  This vignette instantiates

   an object of class \Rclass{CNSet} using a trivial dataset with 3

   files.  As this sample size is too small for estimating copy number

   with the \crlmm{} package, the final section of this vignette loads

   an object created by the analysis of 180 HapMap CEL files

-  (Affymetrix 6.0 platform).  In particular, this object was

-  instantiated by running (1) the \verb+AffymetrixPreprocessCN+

-  vignette and (2) the \verb+copynumber+ vignette.

+  (Affymetrix 6.0 platform).  This object was instantiated by running

+  the \verb+AffyGW+ vignette.

 \end{abstract}

@@ -76,14 +75,14 @@ processing subsets of the markers and/or samples. The functions

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

 values should reduce the RAM required for a particular job.  In

 general, smaller values will increase the 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.  If our dataset

-contained fewer than 500 samples, the \Rfunction{ocSamples} option

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

+a time (when possible) and 200 samples at a time.  If our dataset

+contained fewer than 200 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>>=

-ocProbesets(50e3)

+ocProbesets(100e3)

 ocSamples(200)

 @

@@ -96,11 +95,12 @@ methods:

 \begin{itemize}

-\item[Approach 1:] during the preprocessing of the raw intensities for Illumina and

-  Affymetrix arrays by the the functions \Rfunction{constructInf} and

-  \Rfunction{genotype}, respectively. (The \Rfunction{genotype} calls

-  the non-exported function \Rfunction{constructAffy} to initialize a

-  \Rclass{CNSet} object for Affymetrix platforms.)

+\item[Approach 1:] during the preprocessing of the raw intensities for

+  Illumina and Affymetrix arrays by the the functions

+  \Rfunction{constructInf} and \Rfunction{genotype},

+  respectively. (The \Rfunction{genotype} calls the function

+  \Rfunction{constructAffy} to initialize a \Rclass{CNSet} object for

+  Affymetrix platforms.)

 \item[Approach 2:] by subsetting an existing \Robject{CNSet} object.

   As per usual, the `[' method can be used to extract a subset of

@@ -110,17 +110,18 @@ methods:

 \end{itemize}

 There are important differences in the underlying data representation

-depending on how the object was instantiated.  In particular, objects

-generated by the functions \Rfunction{constructInf} and

-\Rfunction{genotype} store high-dimensional data on disk rather than

-in memory through protocols defined in the \R{} package \ff{}.  For

-instance, the normalized intensities and genotype calls in a

-\Rclass{CNSet}-instance from approach (1) are \ff{}-derived objects.

-By contrast, when such an objected generated by approach (1) is subset

-by the `[' method, an object of the same class is returned but the

-\Rpackage{ff}-derived objects are coerced to ordinary matrices. Note,

-therefore, that both approaches (1) and (2) may involve substantial

-I/O.

+depending on how the \Rclass{CNSet} object was instantiated.  In

+particular, objects generated by the functions

+\Rfunction{constructInf} and \Rfunction{genotype} store

+high-dimensional data on disk rather than in memory through protocols

+defined in the \R{} package \ff{}.  For instance, the normalized

+intensities and genotype calls in a \Rclass{CNSet}-instance from

+approach (1) are \ff{}-derived objects.  By contrast, when such an

+objected generated by approach (1) is subset by the `[' method, an

+object of the same class is returned but the \Rpackage{ff}-derived

+objects are coerced to ordinary matrices. Note, therefore, that both

+approaches (1) and (2) may involve substantial I/O and that (2) should

+be performed judiciously.

 \subsubsection{Approach 1}

@@ -149,7 +150,7 @@ celfiles <- list.celfiles(path, full.names=TRUE)

 Typically, an object of class \Rclass{CNSet} is instantiated as part

 of the preprocessing and genotyping by calling the function

 \Rfunction{genotype}, as illustrated in the

-\verb+AffymtrixPreprocessCN+ vignette.

+\verb+AffyGW+ vignette.

 <<instantiateToyExample,cache=TRUE>>=

 exampleSet <- genotype(celfiles, batch=rep("1", 3), cdfName="genomewidesnp6")

@@ -163,11 +164,11 @@ ldPath()

 @

 One could also instantiate an object of class \Rclass{CNSet} without

-preprocessing/genotyping by calling the non-exported function

+preprocessing/genotyping by calling the exported function

 \Rfunction{constructAffy} directly using the \Rfunction{:::} operator.

 <<constructAffy,eval=FALSE>>=

-tmp <- crlmm:::constructAffy(celfiles, batch=rep("1", 3), cdfName="genomewidesnp6")

+tmp <- constructAffy(celfiles, batch=rep("1", 3), cdfName="genomewidesnp6")

 @

 The \Rfunction{show} method provides a concise summary of the

@@ -300,15 +301,6 @@ platform.  A consequence of keeping the rows of the assay data

 elements the same for all of the statistical summaries is that the

 matrix used to store genotype calls is larger than necessary.

-% only true for affy.  for illumina, we have intensities stored in 'B'

-%Note that NA's are stored in the slot for normalized 'B' allele

-%intensities:

-

-%<<B.NAs>>=

-%np.index <- which(!is.snp)

-%stopifnot(all(is.na(B(cnSet)[np.index, ])))

-%@

-

 \subsubsection{\texttt{batch} and \texttt{batchStatistics}}

 As defined in Leek \textit{et al.} 2010, \textit{Batch effects are

@@ -402,73 +394,10 @@ varLabels(protocolData(exampleSet))

 protocolData(exampleSet)$ScanDate

 @

-

-%\section{Suggested visualizations}

-%

-%\paragraph{SNR.}

-%

-%A histogram of the signal to noise ratio for the HapMap samples:

-%

-%<<plotSnr, fig=TRUE, include=FALSE>>=

-%open(cnSet$SNR)

-%hist(cnSet$SNR[, ], xlab="SNR", main="", breaks=25, col="lightblue", xlim=c(3, max(cnSet$SNR[])))

-%abline(v=5, lty=2, col="grey")

-%text(3,5, label="SNR range for low \n quality arrays", adj=0, col="grey40")

-%@

-%\begin{figure}

-%  \centering

-%  \includegraphics[width=0.6\textwidth]{copynumber-plotSnr}

-%  \caption{Signal to noise ratios for the HapMap samples. SNRs below 5

-%   for the Affymetrix platform are often samples of lower quality.

-%   Such samples will tend to have much more variable estimates of copy

-%   number.}

-%\end{figure}

-

-

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

-%

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

-%versus copy number (vertical axis) for the first sample.  There is

-%less information to estimate copy number at nonpolymorphic loci;

-%improvements to the univariate prediction regions at nonpolymorphic

-%loci are a future area of research. If the \Rpackage{SNPchip} is

-%available, an idiogram can be added to the existing plotting

-%coordinates as indicated in the following example.

-%

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

-%GT.CONF.THR <- 0.8

-%marker.index <- which(chromosome(cnSet) == 1)

-%cn <- totalCopynumber(cnSet, i=marker.index, j=1)

-%x <- position(cnSet)[marker.index]

-%par(las=1, mar=c(4, 5, 4, 2))

-%plot(x, cn, pch=".",

-%     cex=2, xaxt="n", col="grey60", ylim=c(0,6),

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

-%     main=paste(sampleNames(cnSet)[1], ", CHR: 1"))

-%axis(1, at=pretty(x), labels=pretty(x)/1e6)

-%require(SNPchip)

-%invisible(plotCytoband(1, new=FALSE, cytoband.ycoords=c(5.5, 6), label.cytoband=FALSE))

-%@

-

-%\begin{figure}

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

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

-%    chromosome 1 plotted against physical position (x-axis) for one

-%    sample.  Estimates at nonpolymorphic loci are plotted in light

-%    blue.}

-%\end{figure}

-%

-%\clearpage

-%\paragraph{One SNP at a time}

-%

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

-%useful for assessing the biallelic genotype calls.  This section of

-%the vignette is currently under development.

-

 \section{Trouble shooting with a HapMap example}

 This section uses an object of class \Rclass{CNSet} instantiated by

-the \verb+AffymetrixPreprocessCN+ vignette and saved to a local path

+the \verb+AffyGW+ vignette and saved to a local path

 on our computing cluster indicated by the object \Robject{outdir}

 below.  The \verb+copynumber+ vignette was used to fill out the

 \verb+batchStatistics+ slot of the \Robject{cnSet} object.

@@ -490,17 +419,6 @@ invisible(open(cnSet))

 \subsection{Missing values}

-% There are several reasons for estimates of the allele-specific copy

-% number to have missing values (\texttt{NA}'s). This section briefly

-% elaborates on the source of missing values in the HapMap analysis

-% and discusses possible alternatives to reduce the number of missing

-% values.  Note that allele-specific copy number, 'CA' and 'CB', is

-% not saved in the \Robject{cnSet} object.  Rather, the respective

-% accessors calculate 'CA' and 'CB' on the fly from the normalized

-% intensities and from the marker-specific parameter estimates in the

-% linear model.  In general, a missing value arises when the

-% background or slope parameter was not estimated in the linear

-% model.

 Most often, missing values occur when the genotype confidence scores

 for a SNP were below the threshold used by the

 \Robject{crlmmCopynumber} function. For the HapMap analysis, we used a

@@ -129,11 +129,9 @@ directory to store output files and setting the \Rfunction{ldPath}

 function with this path.

 <<path>>=

-if(getRversion() < "2.13.0"){

-	rpath <- getRversion()

-} else rpath <- "trunk"

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

+outdir <- paste("/local_data/r00/crlmm/", getRversion(), "/infrastructure", sep="")

 ldPath(outdir)

+if(!file.exists(outdir)) dir.create(outdir)

 @

@@ -258,6 +256,7 @@ is a SNP can be accessed through accessors defined for the

 \verb+featureData+.

 <<featuredataAccessors>>=

+library(Biobase)

 fvarLabels(exampleSet)

 position(exampleSet)[1:10]

 chromosome(exampleSet)[1:10]

@@ -554,7 +553,7 @@ estimates are missing for a given SNP.

 <<NAchromosomeX>>=

 ## start with first batch

-sample.index <- 1:90

+sample.index <- which(batch(cnSet)==batch(cnSet)[1])

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

 ca.X <- CA(cnSet, i=X.index, j=sample.index)

 missing.caX <- is.na(ca.X)

@@ -22,36 +22,27 @@

 \title{Infrastructure for copy number analysis in \crlmm{}}

 \date{\today}

 \author{Rob Scharpf}

+

 \maketitle

+\begin{abstract}

-\section{Set up}

+  This vignette provides an overview of the \Rclass{CNSet} class and a

+  brief discussion of the underlying infrastructure for large data

+  support with the \Rpackage{ff} package.  This package instantiates

+  an object of class \Rclass{CNSet} using a trivial dataset with 3

+  files.  As this sample size is too small for estimating copy number

+  with the \crlmm{} package, the final section of this vignette loads

+  an object created by the analysis of 180 HapMap CEL files

+  (Affymetrix 6.0 platform).  In particular, this object was

+  instantiated by running (1) the \verb+AffymetrixPreprocessCN+

+  vignette and (2) the \verb+copynumber+ vignette.

-As in the previous vignettes, we load the required libraries and

-specify a path for storing output files.

+\end{abstract}

 <<libraries,results=hide>>=

 library(ff)

 library(crlmm)

-library(cacheSweave)

-require(DNAcopy)

-require(VanillaICE)

-if(getRversion() < "2.13.0"){

-	rpath <- getRversion()

-} else rpath <- "trunk"

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

-@

-

-<<ram>>=

-ldPath(outdir)

-setCacheDir(outdir)

-@

-

-We begin by loading the \Robject{cnSet} object created by the

-\verb+AffymetrixPreprocessCN+ vignette.

-

-<<loadcnset>>=

-if(!exists("cnSet")) load(file.path(outdir, "cnSet.rda"))

 @

 \section{Supported platforms}

@@ -98,31 +89,25 @@ ocSamples(200)

 \section{The \Rclass{CNSet} container}

-The \Rfunction{show} method provides a concise summary of the

-\Robject{cnSet} object.

+\subsection{Instantiating an object of class \Rclass{CNSet}}

-<<show>>=

-cnSet

-@

+An object of class \Rclass{CNSet} can be instantiated by one of two

+methods:

-We briefly outline some of the unique aspects of the

-\Rclass{CNSet}-class using in \crlmm{} that may differ from the more

-standard extensions of the virtual class \Rclass{eSet} defined in the

-\Rpackage{Biobase} package.

+\begin{itemize}

-\paragraph{Instantiating a \Rclass{CNSet}:} An object of class

-\Rclass{CNSet} can be instantiated by one of two methods:

-\begin{enumerate}

-\item during the preprocessing of the raw intensities for Illumina and

+\item[Approach 1:] during the preprocessing of the raw intensities for Illumina and

   Affymetrix arrays by the the functions \Rfunction{constructInf} and

-  \Rfunction{genotype}, respectively.

+  \Rfunction{genotype}, respectively. (The \Rfunction{genotype} calls

+  the non-exported function \Rfunction{constructAffy} to initialize a

+  \Rclass{CNSet} object for Affymetrix platforms.)

-\item by subsetting an existing \Robject{CNSet} object.  As per usual,

-  the `[' method can be used to extract a subset of markers $i$ as in

-  `[i, ]', a subset of samples $j$ as in `[, j]', or a subset of

-  markers $i$ and samples $j$ as in `[i, j]'.

+\item[Approach 2:] by subsetting an existing \Robject{CNSet} object.

+  As per usual, the `[' method can be used to extract a subset of

+  markers $i$ as in `[i, ]', a subset of samples $j$ as in `[, j]', or

+  a subset of markers $i$ and samples $j$ as in `[i, j]'.

-\end{enumerate}

+\end{itemize}

 There are important differences in the underlying data representation

 depending on how the object was instantiated.  In particular, objects

@@ -131,36 +116,106 @@ generated by the functions \Rfunction{constructInf} and

 in memory through protocols defined in the \R{} package \ff{}.  For

 instance, the normalized intensities and genotype calls in a

 \Rclass{CNSet}-instance from approach (1) are \ff{}-derived objects.

-However, when such an object is subset by the `[' method, the data is

-pulled from disk to memory and stored as an ordinary matrix in \R{}.

-To illustrate, the \Robject{cnSet} loaded at the start of this session

-was created by the function \Rfunction{genotype}.  As the data is

-stored on disk rather than in memory, we inspect attributes of the

-object representing the normalized intensities for the A allele by

-first opening the file connection and then extracting the

-\ff{}-derived class of the object as well as where the data is stored

-on disk.

+By contrast, when such an objected generated by approach (1) is subset

+by the `[' method, an object of the same class is returned but the

+\Rpackage{ff}-derived objects are coerced to ordinary matrices. Note,

+therefore, that both approaches (1) and (2) may involve substantial

+I/O.

-<<approach1>>=

-invisible(open(cnSet))

-class(A(cnSet))

-filename(A(cnSet))

+\subsubsection{Approach 1}

+

+To illustrate the first approach, we begin by specifying a local

+directory to store output files and setting the \Rfunction{ldPath}

+function with this path.

+

+<<path>>=

+if(getRversion() < "2.13.0"){

+	rpath <- getRversion()

+} else rpath <- "trunk"

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

+ldPath(outdir)

 @

-Using approach (2), we construct a new \Robject{CNSet} object.

-<<approach2>>=

-cnset.subset <- cnSet[1:5, 1:10]

-invisible(close(cnSet))

-class(A(cnset.subset))

+Next we load the annotation package required, as well as the \R{}

+package \Rpackage{hapmapsnp6} that contains 3 example CEL files.  We

+use the \Rfunction{system.file} function to find the path to the CEL

+files.

+

+<<annotation_metadata>>=

+require(genomewidesnp6Crlmm) & require(hapmapsnp6)

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

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

+@

+

+Typically, an object of class \Rclass{CNSet} is instantiated as part

+of the preprocessing and genotyping by calling the function

+\Rfunction{genotype}, as illustrated in the

+\verb+AffymtrixPreprocessCN+ vignette.

+

+<<instantiateToyExample,cache=TRUE>>=

+exampleSet <- genotype(celfiles, batch=rep("1", 3), cdfName="genomewidesnp6")

 @

-The files with \verb+.ff+ extension should not be removed or

-relocated.  The safest way to move these files if necessary is to

-clone all of the \ff{} objects using the \Rfunction{clone}, followed

-by the \Rfunction{delete} function to remove the original files on

-disk. See the documentation in the \ff{} package for additional

-details.

+Several files with \verb+.ff+ extensions now appear in the directory

+indicated by the \Rfunction{ldPath} function.

+

+<<ldpath>>=

+ldPath()

+@

+

+One could also instantiate an object of class \Rclass{CNSet} without

+preprocessing/genotyping by calling the non-exported function

+\Rfunction{constructAffy} directly using the \Rfunction{:::} operator.

+

+<<constructAffy,eval=FALSE>>=

+tmp <- crlmm:::constructAffy(celfiles, batch=rep("1", 3), cdfName="genomewidesnp6")

+@

+

+The \Rfunction{show} method provides a concise summary of the

+\Robject{exampleSet} object. Note the class of the elements in the

+\verb+batchStatistics+ and \verb+assayData+ slots is indicated in the

+first line of the summary.

+

+<<show>>=

+invisible(open(exampleSet))

+exampleSet

+@

+

+% We briefly outline some of the unique aspects of the

+% \Rclass{CNSet}-class using in \crlmm{} that may differ from the more

+% standard extensions of the virtual class \Rclass{eSet} defined in

+% the \Rpackage{Biobase} package.

+

+As the \verb+assayData+ elements of the \Robject{exampleSet} object

+are stored on disk rather than in memory, we inspect attributes of the

+elements by first opening the file connection using the

+\Rfunction{open}. For example, in the following code we extract the

+normalized intensities for the A allele by first opening the object

+returned by the \Rfunction{A} function.  The name of the file where

+the data is stored on disk is provided by the

+\Rfunction{filename}. Finally, the normalized intensities can be

+pulled from disk to memory by the `[' method.  It can be useful to

+wrap the `[' method by the \Rfunction{as.matrix} to ensure that the

+output is the desired class.

+

+<<approach1>>=

+invisible(open(exampleSet))

+class(A(exampleSet))

+filename(A(exampleSet))

+as.matrix(A(exampleSet)[1:5, ])

+@

+

+\paragraph{Moving \texttt{*.ff} files} Ideally, the files with

+\verb+.ff+ extension should not be moved.  However, if this is not

+possible, the safest way to move these files is to clone all of the

+\ff{} objects using the \Rfunction{clone}, followed by the

+\Rfunction{delete} function to remove the original files on disk. An

+example of the \Rfunction{delete} function is included at the end of

+the \verb+IlluminaPreprocessCN+ vignette. See the documentation for

+the \Rfunction{clone} and \Rfunction{delete} functions in the \ff{}

+package for additional details.

+

 \paragraph{Order of operations:} For \Rclass{CNSet}-instances derived

 by approach (1), users should be mindful of the substantial I/O when

@@ -171,30 +226,49 @@ operation to emphasize the order of operations):

 <<orderOfOperations>>=

 ##inefficient

-A(cnSet[1:5, 1:5])

+##invisible(open(cnSet))

+A(exampleSet[1:5, ])

 ## preferred

-(A(cnSet))[1:5, 1:5]

+(A(exampleSet))[1:5, ]

 @

-\textbf{\texttt{featureData}}

+\subsubsection{Approach 2: using `['}

+

+Here we instantiate a new object of class \Rclass{CNSet} by applying the

+`[' method to an existing object of class \Rclass{CNSet}.

+

+<<approach2>>=

+cnset.subset <- exampleSet[1:5, ]

+@

+

+Note the class of the \verb+batchStatistics+ and \verb+assayData+

+elements of the \Robject{cnset.subset} object printed in the first

+line of the summary.

+

+<<show.subset>>=

+show(cnset.subset)

+@

+

+\subsection{Slots of class \Rclass{CNSet}}

+\subsubsection{\texttt{featureData}}

 Information on physical position, chromosome, and whether the marker

 is a SNP can be accessed through accessors defined for the

 \verb+featureData+.

 <<featuredataAccessors>>=

-fvarLabels(cnSet)

-position(cnSet)[1:10]

-chromosome(cnSet)[1:10]

-is.snp <- isSnp(cnSet)

+fvarLabels(exampleSet)

+position(exampleSet)[1:10]

+chromosome(exampleSet)[1:10]

+is.snp <- isSnp(exampleSet)

 table(is.snp)

 snp.index <- which(is.snp)

 np.index <- which(!is.snp)

-chr1.index <- which(chromosome(cnSet) == 1)

+chr1.index <- which(chromosome(exampleSet) == 1)

 @

-\textbf{\texttt{assayData}}

+\subsubsection{\texttt{assayData}}

 The \verb+assayData+ elements are of the class

 \Rclass{ff\_matrix}/\Rclass{ffdf} or \Rclass{matrix}, depending on how

@@ -203,7 +277,7 @@ the \Rclass{CNSet} object was instantiated.  Elements in the

 function.

 <<assayData>>=

-ls(assayData(cnSet))

+ls(assayData(exampleSet))

 @

 The normalized intensities for the A and B alleles have names

@@ -216,7 +290,7 @@ but can be tranlated to the probability scale using the \R{} function

 \Rfunction{i2p}. For example,

 <<i2p>>=

-scores <- as.matrix(snpCallProbability(cnSet)[1:5, 1:2])

+scores <- as.matrix(snpCallProbability(exampleSet)[1:5, 1:2])

 i2p(scores)

 @

@@ -227,15 +301,16 @@ platform.  A consequence of keeping the rows of the assay data

 elements the same for all of the statistical summaries is that the

 matrix used to store genotype calls is larger than necessary.

-Note that NA's are stored in the slot for normalized 'B' allele

-intensities:

+% only true for affy.  for illumina, we have intensities stored in 'B'

+%Note that NA's are stored in the slot for normalized 'B' allele

+%intensities:

-<<B.NAs>>=

-np.index <- which(!is.snp)

-stopifnot(all(is.na(B(cnSet)[np.index, ])))

-@

+%<<B.NAs>>=

+%np.index <- which(!is.snp)

+%stopifnot(all(is.na(B(cnSet)[np.index, ])))

+%@

-\textbf{\texttt{batch} and \texttt{batchStatistics}}

+\subsubsection{\texttt{batch} and \texttt{batchStatistics}}

 As defined in Leek \textit{et al.} 2010, \textit{Batch effects are

   sub-groups of measurements that have qualitatively different

@@ -260,7 +335,7 @@ to access the \verb+batch+ information on the samples as in the

 following example.

 <<batchAccessor>>=

-table(batch(cnSet))

+batch(exampleSet)

 @

 For the \verb+batchStatistics+ slot, the elements in the environment

@@ -272,7 +347,7 @@ the elements in the environment can be list using the \R{} function

 \Rfunction{ls}.

 <<batchStatistics>>=

-ls(batchStatistics(cnSet))

+ls(batchStatistics(exampleSet))

 @

 Currently, the batch-specific summaries are stored to allow some

@@ -292,7 +367,7 @@ intended to be accessed directly by the user.

 %median variance (across markers) within a batch. ) The elements of the

 %slot can be listed as follows.

-\textbf{\texttt{phenoData}}

+\subsubsection{\texttt{phenoData}}

 Sample-level summaries obtained during the preprocessing/genotyping

 steps include skew (SKW), the signal to noise ratio (SNR), and gender

@@ -303,26 +378,29 @@ approach (1), these elements are of the class

 \Rclass{ff\_vector}.

 <<phenodataAccessors>>=

-varLabels(cnSet)

-class(cnSet$gender)

-invisible(open(cnSet$gender))

-cnSet$gender

+varLabels(exampleSet)

+class(exampleSet$gender)

+invisible(open(exampleSet$gender))

+exampleSet$gender

 @

 The `[' methods without arguments can be used to coerce to a vector.

 <<vector>>=

-cnSet$gender[]

-invisible(close(cnSet$gender))

+c("male", "female")[exampleSet$gender[]]

+invisible(close(exampleSet$gender))

 @

-\textbf{\texttt{protocolData}}

+

+

+

+\subsubsection{\texttt{protocolData}}

 The scan date of the arrays are stored in the \verb+protocolData+.

 <<protocolData>>=

-varLabels(protocolData(cnSet))

-protocolData(cnSet)$ScanDate[1:10]

+varLabels(protocolData(exampleSet))

+protocolData(exampleSet)$ScanDate

 @

@@ -388,7 +466,28 @@ protocolData(cnSet)$ScanDate[1:10]

 %useful for assessing the biallelic genotype calls.  This section of

 %the vignette is currently under development.

-\section{Trouble shooting}

+\section{Trouble shooting with a HapMap example}

+

+This section uses an object of class \Rclass{CNSet} instantiated by

+the \verb+AffymetrixPreprocessCN+ vignette and saved to a local path

+on our computing cluster indicated by the object \Robject{outdir}

+below.  The \verb+copynumber+ vignette was used to fill out the

+\verb+batchStatistics+ slot of the \Robject{cnSet} object.

+

+<<path>>=

+if(getRversion() < "2.13.0"){

+	rpath <- getRversion()

+} else rpath <- "trunk"

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

+@

+

+Next, we load the \Robject{cnSet} object.

+

+<<loadcnset>>=

+if(!exists("cnSet")) load(file.path(outdir, "cnSet.rda"))

+invisible(open(cnSet))

+@

+

 \subsection{Missing values}

@@ -406,12 +505,15 @@ protocolData(cnSet)$ScanDate[1:10]

 Most often, missing values occur when the genotype confidence scores

 for a SNP were below the threshold used by the

 \Robject{crlmmCopynumber} function. For the HapMap analysis, we used a

-confidence threshold of 0.80 (the default).

+confidence threshold of 0.80 (the default).  In the following code, we

+assess NA's appearing for the raw copy number estimates for the first

+10 samples.

 <<NA>>=

 GT.CONF.THR <- 0.80

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

-sample.index <- which(batch(cnSet) == "C")

+sample.index <- 1:10

+ct <- totalCopynumber(cnSet, i=autosome.index, j=sample.index)

 ca <- CA(cnSet, i=autosome.index, j=sample.index)

 cb <- CB(cnSet, i=autosome.index, j=sample.index)

 missing.ca <- is.na(ca)

@@ -427,9 +529,9 @@ threshold specified in \Robject{crlmmCopynumber}.

 <<NAconfidenceScores>>=

 if(nmissing.ca > 0){

-	invisible(open(snpCallProbability(cnSet)))

+	##invisible(open(snpCallProbability(cnSet)))

 	gt.conf <- i2p(snpCallProbability(cnSet)[autosome.index, sample.index])

-	invisible(close(snpCallProbability(cnSet)))

+	##invisible(close(snpCallProbability(cnSet)))

 	below.thr <- gt.conf < GT.CONF.THR

 	index.allbelow <- as.integer(which(rowSums(below.thr) == length(sample.index)))

 	nmissingBecauseOfGtThr <- length(index.allbelow) * length(sample.index)

@@ -451,13 +553,15 @@ missing values to the number of samples to check whether all of the

 estimates are missing for a given SNP.

 <<NAchromosomeX>>=

+## start with first batch

+sample.index <- 1:90

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

 ca.X <- CA(cnSet, i=X.index, j=sample.index)

 missing.caX <- is.na(ca.X)

 (nmissing.caX <- sum(missing.caX))

 missing.snp.index <- which(rowSums(missing.caX) == length(sample.index))

 index <- which(rowSums(missing.caX) == length(sample.index))

-length(index)*length(sample.index)/nmissing.caX

+##length(index)*length(sample.index)/nmissing.caX

 @

 From the above codechunk, we see that \Sexpr{length(index)} SNPs have

@@ -468,21 +572,15 @@ from SNPs in which either the men or the women had confidence scores

 that were all below the threshold.

 <<NAcrlmmConfidenceScore>>=

-if(nmissing.caX > 0){

-	invisible(open(snpCallProbability(cnSet)))

-	gt.conf <- i2p(snpCallProbability(cnSet)[X.index, sample.index])

-	invisible(close(snpCallProbability(cnSet)))

-	below.thr <- gt.conf < GT.CONF.THR

-	F <- which(cnSet$gender[sample.index] == 2)

-	M <- which(cnSet$gender[sample.index] == 1)

-	##nus <- nuA(cnSet)[X.index, "C"]

-	index.allbelowF <- as.integer(which(rowSums(below.thr[, F]) == length(F)))

-	index.allbelowM <- as.integer(which(rowSums(below.thr[, M]) == length(M)))

-	index.allbelow <- as.integer(unique(c(index.allbelowF, index.allbelowM)))

-	##all.equal(index, index.allbelow)

-	## or calculate the proportion of missing effected by low crlmm confidence

-	sum(index.allbelow %in% index)/length(index)

-}

+invisible(open(cnSet$gender))

+F <- which(cnSet$gender[sample.index] == 2)

+M <- which(cnSet$gender[sample.index] == 1)

+gt.conf <- i2p(snpCallProbability(cnSet)[X.index, sample.index])

+below.thr <- gt.conf < GT.CONF.THR

+index.allbelowF <- as.integer(which(rowSums(below.thr[, F]) == length(F)))

+index.allbelowM <- as.integer(which(rowSums(below.thr[, M]) == length(M)))

+all.equal(index.allbelowF, index.allbelowM)

+all.equal(as.integer(index.allbelowF), as.integer(index))

 @

 For nonpolymorphic loci, the genotype confidence scores are irrelevant

@@ -519,8 +617,9 @@ detection and removal may be explored in the future.

 Copy number estimates for other chromosomes, such as mitochondrial and

 chromosome Y, are not currently available in \crlmm{}.

-<<echo=FALSE>>=

+<<close>>=

 invisible(close(cnSet))

+invisible(close(cnSet$gender))

 @

@@ -1,4 +1,4 @@

-%\VignetteIndexEntry{Infrastructure for crlmm copy number Vignette}

+%\VignetteIndexEntry{Infrastructure for copy number analysis}

 %\VignetteDepends{crlmm, genomewidesnp6Crlmm}

 %\VignetteKeywords{crlmm, SNP 6}

 %\VignettePackage{crlmm}

new file mode 100644

@@ -0,0 +1,540 @@

+%\VignetteIndexEntry{Infrastructure for crlmm copy number Vignette}

+%\VignetteDepends{crlmm, genomewidesnp6Crlmm}

+%\VignetteKeywords{crlmm, SNP 6}

+%\VignettePackage{crlmm}

+\documentclass{article}

+\usepackage{graphicx}

+\usepackage{natbib}

+\usepackage{url}

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

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

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

+\newcommand{\Robject}[1]{{\texttt{#1}}}

+\newcommand{\Rpackage}[1]{{\textsf{#1}}}

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

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

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

+\newcommand{\crlmm}{\Rpackage{crlmm}}

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

+\usepackage[margin=1in]{geometry}

+

+\begin{document}

+\title{Infrastructure for copy number analysis in \crlmm{}}

+\date{\today}

+\author{Rob Scharpf}

+\maketitle

+

+

+\section{Set up}

+

+As in the previous vignettes, we load the required libraries and

+specify a path for storing output files.

+

+<<libraries,results=hide>>=

+library(ff)

+library(crlmm)

+library(cacheSweave)

+require(DNAcopy)

+require(VanillaICE)

+if(getRversion() < "2.13.0"){

+	rpath <- getRversion()

+} else rpath <- "trunk"

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

+@

+

+<<ram>>=

+ldPath(outdir)

+setCacheDir(outdir)

+@

+

+We begin by loading the \Robject{cnSet} object created by the

+\verb+AffymetrixPreprocessCN+ vignette.

+

+<<loadcnset>>=

+if(!exists("cnSet")) load(file.path(outdir, "cnSet.rda"))

+@

+

+\section{Supported platforms}

+

+The supported Affymetrix and 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. For

+Affymetrix, the 5.0 and 6.0 platforms are supported and the

+corresponding annotation packages are \Rpackage{genomewidesnp5Crlmm}

+and \Rpackage{genomewidesnp6Crlmm}.  Supported Infinium platforms are

+listed in the following code chunk.

+

+<<supportedPlatforms>>=

+pkgs <- annotationPackages()

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

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

+@

+

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

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

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

+writing data to disk instead of keeping data in memory.  As the

+functions for preprocessing, genotyping, and copy number estimation do

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

+general, smaller values will increase the 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.  If our dataset

+contained fewer than 500 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>>=

+ocProbesets(50e3)

+ocSamples(200)

+@

+

+\section{The \Rclass{CNSet} container}

+

+The \Rfunction{show} method provides a concise summary of the

+\Robject{cnSet} object.

+

+<<show>>=

+cnSet

+@

+

+We briefly outline some of the unique aspects of the

+\Rclass{CNSet}-class using in \crlmm{} that may differ from the more

+standard extensions of the virtual class \Rclass{eSet} defined in the

+\Rpackage{Biobase} package.

+

+\paragraph{Instantiating a \Rclass{CNSet}:} An object of class