@@ -200,8 +200,7 @@ Estimation of log R ratios and B allele frequencies from an object of class \Rcl

 crlmmCopynumber(cnSet, fit.linearModel=FALSE)

 @

-<<oligoSnpSet,eval=FALSE>>=

-library(VanillaICE)

+<<oligoSnpSet>>=

 open(cnSet)

 oligoSetList <- BafLrrSetList(cnSet)

 close(cnSet)

@@ -214,7 +213,7 @@ oligoSetList[[1]]

 \noindent Log R ratios and B allele frequences can be retrieved by the

 accessors \Rfunction{lrr} and \Rfunction{baf}, respectively.

-<<testEqual,eval=FALSE>>=

+<<testEqual>>=

 lrrList <- lrr(oligoSetList)

 class(lrrList)

 dim(lrrList[[1]]) ## log R ratios for chromosome 1.

@@ -225,4 +225,3 @@ dim(bafList[[1]]) ## B allele frequencies for chromosome 1

 If the \Rfunction{crlmmCopynumber} function is not run prior to the

 \Rfunction{BafLrrSetList} construction, the log R ratios and BAFs will

 be initialized as NAs.

-

@@ -53,15 +53,24 @@ print(histogram(~snr,

 \section{Copy number estimation}

-As described in \cite{Scharpf2011}, the CRLMM-CopyNumber algorithm

-fits a linear model to the normalized intensities stratified by the

-diallic genotype call.  The intercept and slope from the linear model

-are both SNP- and batch-specific.  The implementation in the \crlmm{}

-package is encapsulated by the function \Rfunction{crlmmCopynumber}

-that, using the default settings, can be called by passing a single

-object of class \Rclass{CNSet}.  See the appropriate

-preprocessing/genotyping vignette for the construction of an object of

-class \Rclass{CNSet}.

+There are two ways to obtain marker-level estimates of copy number

+that are supported by \Rpackage{crlmm}. One approach is to fit a

+linear model to the normalized intensities stratified by the diallelic

+genotype call at each SNP, as described in \cite{Scharpf2011}.

+Another alternative is to compute log R ratio and B allele

+frequencies. The latter is often better supported by downstream hidden

+Markov models such as those in the \Rpackage{VanillaICE} package. We

+describe each approach in the following two sections.

+

+\subsection{Linear model for normalized intensities}

+

+In this section, we fit a linear model to the normalized intensities

+stratified by the diallic genotype call.  The intercept and slope from

+the linear model are both SNP- and batch-specific.  The implementation

+in the \crlmm{} package is encapsulated by the function

+\Rfunction{crlmmCopynumber} that, using the default settings, can be

+called by passing a single object of class \Rclass{CNSet}.

+

 <<LDS_copynumber,cache=TRUE>>=

 crlmmCopynumber(cnSet)

 @

@@ -105,12 +114,7 @@ described in the \R{} package \Rpackage{ff}. The value returned by

 \Rfunction{crlmmCopynumber} is \verb+TRUE+, indicating that the files

 on disk have been successfully updated.  Note that while the

 \Robject{cnSet} object is unchanged, the values on disk are different.

-On the other hand, subsetting the \Robject{cnSet} with the `[' method

-coerces all of the elements to class \Rclass{matrix}. The

-batch-specific summaries are now ordinary matrices stored in RAM. The

-object returned by \Robject{crlmmCopynumber} is an object of class

-\Rclass{CNSet} with the matrices in the \verb+batchStatistics+ slot

-updated.

+

 <<chr1index>>=

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

@@ -127,20 +131,10 @@ for(i in seq_along(nms)) cls[i] <- class(batchStatistics(cnSet2)[[nms[i]]])[1]

 all(cls == "matrix")

 @

-<<matrix_copynumber,cache=TRUE>>=

-cnSet3 <- crlmmCopynumber(cnSet2)

-class(cnSet3)

-@

-

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

 rm(cnSet2); gc()

 @

-\subsection{Marker-specific estimates}

-%\subsection{Raw copy number}

-

-%\paragraph{Log R ratios and B allele frequencies.}

-

 \paragraph{Raw total copy number.}

 Several functions are available that will compute relatively quickly

 the allele-specific, \emph{raw} copy number estimates. At allele $k$,

@@ -186,38 +180,49 @@ cb <- CB(cnSet, i=snp.index, j=1:2)

 \subsection{Container for log R ratios and B allele frequencies}

+\label{sec:lrrBaf}

 A useful container for storing the \crlmm{} genotypes, genotype

 confidence scores, and the total or relative copy number at each

-marker is the \Rclass{oligoSetList} class.  Coercion of a

-\Rclass{CNSet} object to a \Rfunction{oligoSnpSet} object can be

-acheived by using the function \Rfunction{constructOligoSetFrom} as

+marker is the \Rclass{BafLrrSetList} class.  Coercion of a

+\Rclass{CNSet} object to a \Rfunction{BafLrrSetList} object can be

+acheived by the function \Rfunction{BafLrrSetList} as

 illustrated below. Users should note that if the \verb+assayData+

 elements in the \Rclass{CNSet} instance are \Rpackage{ff} objects, the

-\verb+assayData+ elements of each element in the \Rclass{oligoSetList}

+\verb+assayData+ elements of each element in the \Rclass{BafLrrSetList}

 object will be \Rpackage{ff}-dervied objects (a new

 \verb+total_cn*.ff+ file will be created in the \verb+ldPath()+

-directory).

+directory). The following code-chunk is not evalutated.

+

+Estimation of log R ratios and B allele frequencies from an object of class \Rclass{CNSet} instantiated from genotyping Affymetrix CEL files or Illumina iDat files requires running the \R{} function \Rfunction{crlmmCopynumber} as a preliminary step.  Specifying \texttt{fit.linearModel=FALSE} will speed up computation as this function will by default compute summary statistics unnecessary for computing BAFs and log R ratios.

-<<oligoSnpSet>>=

+<<callingCrlmmCopynumber,eval=FALSE>>=

+crlmmCopynumber(cnSet, fit.linearModel=FALSE)

+@

+

+<<oligoSnpSet,eval=FALSE>>=

 library(VanillaICE)

-open(cnSet3)

-oligoSetList <- constructOligoSetListFrom(cnSet3, batch.name=batch(cnSet3)[1])

-close(cnSet3)

+open(cnSet)

+oligoSetList <- BafLrrSetList(cnSet)

+close(cnSet)

 show(oligoSetList)

 class(oligoSetList)

 ## oligoSnpSet of first chromosome

 oligoSetList[[1]]

 @

-\noindent Note that log R ratios stored in the \Robject{oligoSnpSet}

-object can be retrieved by the \Rfunction{copyNumber} accessor.  B

-allele frequences are retrieved by the \Rfunction{baf} accessor.

+\noindent Log R ratios and B allele frequences can be retrieved by the

+accessors \Rfunction{lrr} and \Rfunction{baf}, respectively.

-<<testEqual>>=

-lrrList <- copyNumber(oligoSetList)

+<<testEqual,eval=FALSE>>=

+lrrList <- lrr(oligoSetList)

 class(lrrList)

 dim(lrrList[[1]]) ## log R ratios for chromosome 1.

 bafList <- baf(oligoSetList)

 dim(bafList[[1]]) ## B allele frequencies for chromosome 1

 @

+

+If the \Rfunction{crlmmCopynumber} function is not run prior to the

+\Rfunction{BafLrrSetList} construction, the log R ratios and BAFs will

+be initialized as NAs.

+

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

   Affymetrix 5.0 and 6.0 platforms, as well as several Illumina

   platforms.  This vignette assumes that the arrays have already been

   successfully preprocessed and genotyped as per the instructions in

-  the \verb+AffymetrixPreprocessCN+ and \verb+IlluminaPreprocessCN+

+  the \verb+AffyGW+ and \verb+IlluminaPreprocessCN+

   vignettes for the Affymetrix and Illumina platforms,

   respectively. While this vignette uses Affymetrix 6.0 arrays for

   illustration, the steps at this point are identical for both

@@ -63,7 +63,7 @@ object of class \Rclass{CNSet}.  See the appropriate

 preprocessing/genotyping vignette for the construction of an object of

 class \Rclass{CNSet}.

 <<LDS_copynumber,cache=TRUE>>=

-(cnSet.updated <- crlmmCopynumber(cnSet))

+crlmmCopynumber(cnSet)

 @

 The following steps were performed by the \Rfunction{crlmmCopynumber}

@@ -154,7 +154,7 @@ formally, \newcommand{\A}{A} \newcommand{\B}{B}

     \phi}_{k,ip}}\left(I_{k,ijp}-{\hat \nu}_{k,ip}\right), ~0\right\}

 \mbox{~for~} k \in \{\A, \B\}.

 \end{eqnarray}

-\noindent See \cite{Scharpf2010} for details.

+\noindent See \cite{Scharpf2011} for details.

 The function \Rfunction{totalCopynumber} translates the normalized

 intensities to an estimate of raw copy number by adding the

@@ -175,42 +175,49 @@ dim(tmp2)

 Alternatively, the functions \Rfunction{CA} and \Rfunction{CB} compute

 the allele-specific copy number.  For instance, the following code

-chunk computes the allele-specific summaries at all polymorphic loci.

+chunk computes the allele-specific summaries at all polymorphic loci

+for the first 2 samples.

 <<ca>>=

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

-ca <- CA(cnSet, i=snp.index, j=seq_len(ncol(cnSet)))

-cb <- CB(cnSet, i=snp.index, j=seq_len(ncol(cnSet)))

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

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

 @

-\subsection{Container for raw copy number}

+\subsection{Container for log R ratios and B allele frequencies}

 A useful container for storing the \crlmm{} genotypes, genotype

-confidence scores, and the total copy number at each marker is the

-\Rclass{oligoSnpSet} class.  Coercion of a \Rclass{CNSet} object to a

-\Rfunction{oligoSnpSet} object can be acheived by using the method

-\Rfunction{as} (as illustrated below). Users should note that if the

-\verb+assayData+ elements in the \Rclass{CNSet} instance are

-\Rpackage{ff} objects, the \verb+assayData+ elements of the

-instantiated \Rclass{oligoSnpSet} will also be \Rpackage{ff}-dervied

-objects (a new \verb+total_cn*.ff+ file will be created in the

-\verb+ldPath()+ directory).

+confidence scores, and the total or relative copy number at each

+marker is the \Rclass{oligoSetList} class.  Coercion of a

+\Rclass{CNSet} object to a \Rfunction{oligoSnpSet} object can be

+acheived by using the function \Rfunction{constructOligoSetFrom} as

+illustrated below. Users should note that if the \verb+assayData+

+elements in the \Rclass{CNSet} instance are \Rpackage{ff} objects, the

+\verb+assayData+ elements of each element in the \Rclass{oligoSetList}

+object will be \Rpackage{ff}-dervied objects (a new

+\verb+total_cn*.ff+ file will be created in the \verb+ldPath()+

+directory).

 <<oligoSnpSet>>=

+library(VanillaICE)

 open(cnSet3)

-oligoSet <- as(cnSet3, "oligoSnpSet")

+oligoSetList <- constructOligoSetListFrom(cnSet3, batch.name=batch(cnSet3)[1])

 close(cnSet3)

-class(copyNumber(oligoSet))

+show(oligoSetList)

+class(oligoSetList)

+## oligoSnpSet of first chromosome

+oligoSetList[[1]]

 @

-\noindent Note that the raw copy number estimates stored in the

-\Robject{oligoSnpSet} object can be retrieved by the

-\Rfunction{copyNumber} accessor and is equivalent to that returned by

-the \Rfunction{totalCopynumber} function defined over the same row and

-column indices.

+\noindent Note that log R ratios stored in the \Robject{oligoSnpSet}

+object can be retrieved by the \Rfunction{copyNumber} accessor.  B

+allele frequences are retrieved by the \Rfunction{baf} accessor.

 <<testEqual>>=

-total.cn3 <- totalCopynumber(cnSet3, i=1:nrow(cnSet3), j=seq_len(ncol(cnSet3)))

-all.equal(copyNumber(oligoSet), total.cn3)

-@

\ No newline at end of file

+lrrList <- copyNumber(oligoSetList)

+class(lrrList)

+dim(lrrList[[1]]) ## log R ratios for chromosome 1.

+bafList <- baf(oligoSetList)

+dim(bafList[[1]]) ## B allele frequencies for chromosome 1

+@

@@ -1,35 +1,3 @@

-%\VignetteIndexEntry{crlmm copy number Vignette}

-%\VignetteDepends{crlmm, genomewidesnp6Crlmm}

-%\VignetteKeywords{crlmm, SNP 6}

-%\VignettePackage{crlmm}

-\documentclass{article}

-\usepackage{graphicx}

-\usepackage{natbib}

-\usepackage{url}

-\usepackage{amsmath,amssymb}

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

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

-

-\begin{document}

-\title{Copy number estimation with \Rpackage{crlmm}}

-\date{\today}

-\author{Rob Scharpf}

-\maketitle

-

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

-options(continue=" ", width=70)

-@

-

-\begin{abstract}

-

   Copy number routines in the \crlmm{} package are available for

   Affymetrix 5.0 and 6.0 platforms, as well as several Illumina

   platforms.  This vignette assumes that the arrays have already been

@@ -44,35 +12,6 @@ options(continue=" ", width=70)

   \crlmm{} package is available from the author's website:

   \url{http://www.biostat.jhsph.edu/~rscharpf/crlmmCompendium/index.html}.

-\end{abstract}

-

-\section{Set up}

-

-<<libraries,results=hide>>=

-library(ff)

-library(crlmm)

-library(lattice)

-library(cacheSweave)

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

-ocProbesets(150e3)

-ocSamples(200)

-@

-

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

-\verb+AffymetrixPreprocessCN+ vignette.

-

-<<loadcnset>>=

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

-@

-

 \textbf{Limitations:} While a minimum number of samples is not

 required for preprocessing and genotyping, copy number estimation in

@@ -100,24 +39,16 @@ quality.  The following code chunk makes a histogram of the SNR values

 for the HapMap samples.

 <<snr,fig=TRUE,include=FALSE,width=6, height=4>>=

+library(lattice)

 invisible(open(cnSet$SNR))

 snr <- cnSet$SNR[]

 close(cnSet$SNR)

-print(histogram(~snr, panel=function(...){

-	panel.histogram(...)

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

-},

-		breaks=25,xlim=c(4.5,10), xlab="SNR"))

+print(histogram(~snr,

+		panel=function(...){

+			panel.histogram(...)},

+		breaks=25, xlim=range(snr), xlab="SNR"))

 @

-\begin{figure}[t!]

-  \begin{center}

-  \includegraphics[width=0.6\textwidth]{copynumber-snr.pdf}

-  \caption{The signal to noise ratio (SNR) for 180 HapMap samples. For

-    Affymetrix platforms, SNR values below 5 can indicate possible

-    problems with sample quality.  In some circumstances, it may be

-    more helpful to exclude samples with poor DNA quality.}

-\end{center}

-\end{figure}

+

 \section{Copy number estimation}

@@ -149,29 +80,17 @@ computed for subsets of the markers to reduce the required RAM. Note

 that the value returned by the \Rfunction{crlmmCopynumber} function in

 the above example is \verb+TRUE+.  The reason the function returns

 \verb+TRUE+ in the above example is that the elements of the

-\verb+batchStatistics+ slot have the class \Rclass{ff\_matrix}.  Rather

-than keep the statistical summaries in memory, the summaries are

-written to files on disk using protocols described in the

+\verb+batchStatistics+ slot have the class \Rclass{ff\_matrix}.

+Rather than keep the statistical summaries in memory, the summaries

+are written to files on disk using protocols described in the

 \Rpackage{ff} package. Hence, while the \Robject{cnSet} object itself

 is unchanged as a result of the \Rfunction{crlmmCopynumber} function,

 the data on disk is updated accordingly.  Users that are interested in

 accessing these low-level summaries can refer to the

-\verb+Infrastructure+ vignette.  Computation of the raw copy number

-estimates for each allele is described in the following section.

-

-For users that are interested in the analysis of a specific chromosome

-(subset of markers) or a s

-

-pointers to files on disk, are stored in the \verb+batchStatistics+

-slot of the class \Rclass{CNSet}.  Using the default settings for the

-\Rfunction{crlmmCopynumber} function, only an object of class

-\Rclass{CNSet} is required.  %(See \verb+AffyPreprocessCN+ or

-%\verb+IlluminaPreprocessCN+ vignettes for details.)

-

-Note that depends on whether the elements of the \verb+batchStatistics+

-slot are \Robject{ff} objects or ordinary matrices.  In this example,

-the elements of \verb+batchStatistics+ have the class

-\Rclass{ff\_matrix}.

+\verb+Infrastructure+ vignette.  Note that the data structure depends

+on whether the elements of the \verb+batchStatistics+ slot are

+\Robject{ff} objects or ordinary matrices.  In this example, the

+elements of \verb+batchStatistics+ have the class \Rclass{ff\_matrix}.

 <<classes>>=

 nms <- ls(batchStatistics(cnSet))

@@ -186,10 +105,8 @@ described in the \R{} package \Rpackage{ff}. The value returned by

 \Rfunction{crlmmCopynumber} is \verb+TRUE+, indicating that the files

 on disk have been successfully updated.  Note that while the

 \Robject{cnSet} object is unchanged, the values on disk are different.

-

-

-\noindent On the other hand, subsetting the \Robject{cnSet} with the

-`[' method coerces all of the elements to class \Rclass{matrix}. The

+On the other hand, subsetting the \Robject{cnSet} with the `[' method

+coerces all of the elements to class \Rclass{matrix}. The

 batch-specific summaries are now ordinary matrices stored in RAM. The

 object returned by \Robject{crlmmCopynumber} is an object of class

 \Rclass{CNSet} with the matrices in the \verb+batchStatistics+ slot

@@ -219,14 +136,18 @@ class(cnSet3)

 rm(cnSet2); gc()

 @

-\subsection{Raw copy number}

+\subsection{Marker-specific estimates}

+%\subsection{Raw copy number}

+

+%\paragraph{Log R ratios and B allele frequencies.}

+\paragraph{Raw total copy number.}

 Several functions are available that will compute relatively quickly

 the allele-specific, \emph{raw} copy number estimates. At allele $k$,

 marker $i$, sample $j$, and batch $p$, the estimate of allele-specific

 copy number is computed by subtracting the estimated background from

-the normalized intensity and scaling by the slope coefficient. More formally,

-\newcommand{\A}{A} \newcommand{\B}{B}

+the normalized intensity and scaling by the slope coefficient. More

+formally, \newcommand{\A}{A} \newcommand{\B}{B}

 \begin{eqnarray}

   \label{eq:cnK}

 {\hat c}_{k,ijp} &=& \mbox{max}\left\{\frac{1}{{\hat

@@ -246,9 +167,9 @@ at all markers for the first 2 samples, and the total copy number for

 chromosome 20 for the first 50 samples.

 <<totalCopynumber>>=

-tmp <- totalCopynumber(cnSet, i=1:nrow(cnSet), j=1:2)

+tmp <- totalCopynumber(cnSet, i=seq_len(nrow(cnSet)), j=1:2)

 dim(tmp)

-tmp2 <- totalCopynumber(cnSet, i=which(chromosome(cnSet) == 20), j=1:50)

+tmp2 <- totalCopynumber(cnSet, i=which(chromosome(cnSet) == 20), j=seq_len(ncol(cnSet)))

 dim(tmp2)

 @

@@ -258,98 +179,12 @@ chunk computes the allele-specific summaries at all polymorphic loci.

 <<ca>>=

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

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

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

-@

-

-\noindent Note the equivalence of the following calculations.

-

-<<totalCopynumber>>=

-ct <- ca+cb

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

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

-@

-

-At nonpolymorphic loci, \Rfunction{CA} function returns the total copy

-number and, by construction, the \Rfunction{CB} function returns 0.

-

-<<nonpolymorphicAutosomes>>=

-marker.index <- which(!isSnp(cnSet))

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

-@

-

-In the following code chunk, we extract estimates of the total copy

-number at nonpolymorphic markers on chromosome X.

-

-<<npx>>=

-set.seed(123)

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

-@

-

-\noindent Again, the function \Rfunction{totalCopynumber} is

-equivalent.

-

-<<npx2>>=

-cnX <- cbind(cn.M, cn.F)

-cnX2 <- totalCopynumber(cnSet, i=npx.index, j=c(M, F))

-stopifnot(all.equal(cnX, cnX2))

+ca <- CA(cnSet, i=snp.index, j=seq_len(ncol(cnSet)))

+cb <- CB(cnSet, i=snp.index, j=seq_len(ncol(cnSet)))

 @

-<<nonpolymorphicX,fig=TRUE,include=FALSE, echo=FALSE, eval=FALSE>>=

-df <- data.frame(cn=as.numeric(cnX), id=factor(rep(sampleNames(cnSet)[c(M,F)], each=nrow(cnX)), levels=sampleNames(cnSet)[c(M,F)], ordered=T))

-library(RColorBrewer)

-cols <- brewer.pal(8, "Accent")[c(1, 5)]

-print(bwplot(cn~id,df, panel=function(x,y,...){

-	panel.grid(v=-10,h=0)

-	panel.bwplot(x,y, ...)

-	panel.abline(h=1:2, col="grey70",lwd=2)

-},

-	scales=list(x=list(rot=90)), cex=0.5, ylab="total copy number", main="Chr X SNPs",

-       fill=cols[cnSet$gender[c(M,F)]], key=mykey))

-@

-

-%\begin{figure}[t!]

-% \centering

-% \includegraphics[width=0.5\textwidth]{copynumber-nonpolymorphicX}

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

-%   X (5 men, 5 women).  We assume that the median copy number across

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

-%\end{figure}

-

-Polymorphic markers on chromosome X:

-<<polymorphicX,fig=TRUE,include=FALSE>>=

-library(RColorBrewer)

-cols <- brewer.pal(8, "Accent")[c(1, 5)]

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

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

-df <- data.frame(cn=as.numeric(cnX), id=factor(rep(sampleNames(cnSet)[c(M,F)], each=length(X.markers)), levels=sampleNames(cnSet)[c(M,F)], ordered=T))

-mykey <- simpleKey(c("male", "female"), points=FALSE,col=cols)

-print(bwplot(cn~id,df, panel=function(x,y,...){

-	panel.grid(v=-10,h=0)

-	panel.bwplot(x,y, ...)

-	panel.abline(h=1:2, col="grey70",lwd=2)

-},

-	scales=list(x=list(rot=90)), cex=0.5, ylab="total copy number", main="Chr X SNPs",

-       fill=cols[cnSet$gender[c(M,F)]], key=mykey))

-@

-\begin{figure}[t!]

- \centering

- \includegraphics[width=0.5\textwidth]{copynumber-polymorphicX}

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

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

-

-\subsection{A container for raw copy number}

+\subsection{Container for raw copy number}

 A useful container for storing the \crlmm{} genotypes, genotype

 confidence scores, and the total copy number at each marker is the

@@ -362,28 +197,6 @@ instantiated \Rclass{oligoSnpSet} will also be \Rpackage{ff}-dervied

 objects (a new \verb+total_cn*.ff+ file will be created in the

 \verb+ldPath()+ directory).

-%The advantage of this approach is that the raw copy number estimates

-%can be computed for all markers and all samples simultaneously without

-%excessive RAM.  The disadvantage is that such a step creates

-%additional I/O for reading and writing the raw copy number estimates

-%from disk.  For very large data sets (1000+ samples), the time

-%required for I/O may outweigh the benefits.  The alternative is to

-%perform downstream processing (such as smoothing) on the total copy

-%number for subsets of markers and/or samples.  As indicated

-%previously, the subset `[' method applied to an object of class

-%\Robject{CNSet} will return an object of the same class with simple

-%matrices instead of \Rpackage{ff} objects.  Therefore, one can

-%instantiate a \Rpackage{oligoSnpSet} object from a \Robject{CNSet}

-%object without the I/O overhead for storing the total copy number

-%estimates by first subsetting the \Rclass{CNSet} object.  One could

-%then smooth the raw copy number estimates and store a summary copy

-%number for the interval, thereby substantially reducing the

-%dimensionality of the copy number estimates.  Examples of smoothing

-%applications such as circular binary segmentation implemented in the

-%\Rpackage{DNAcopy} package or a hidden Markov model implemented in the

-%\Rpackage{VanillaICE} package are provided in the \verb+SmoothingCN+

-%vignette.

-

 <<oligoSnpSet>>=

 open(cnSet3)

 oligoSet <- as(cnSet3, "oligoSnpSet")

@@ -398,19 +211,6 @@ the \Rfunction{totalCopynumber} function defined over the same row and

 column indices.

 <<testEqual>>=

-total.cn3 <- totalCopynumber(cnSet3, i=1:nrow(cnSet3), j=1:ncol(cnSet3))

+total.cn3 <- totalCopynumber(cnSet3, i=1:nrow(cnSet3), j=seq_len(ncol(cnSet3)))

 all.equal(copyNumber(oligoSet), total.cn3)

-@

-

-

-\section{Session information}

-<<sessionInfo, results=tex>>=

-toLatex(sessionInfo())

-@

-

-%\begin{bibliography}

-  \bibliographystyle{plain}

-  \bibliography{refs}

-%\end{bibliography}

-

-\end{document}

+@

\ No newline at end of file

@@ -38,7 +38,7 @@ options(continue=" ", width=70)

   vignettes for the Affymetrix and Illumina platforms,

   respectively. While this vignette uses Affymetrix 6.0 arrays for

   illustration, the steps at this point are identical for both

-  platforms.  See \citep{Scharpf2010} for details regarding the

+  platforms.  See \citep{Scharpf2011} for details regarding the

   methodology implemented in \crlmm{} for copy number analysis.  In

   addition, a compendium describing copy number analysis using the

   \crlmm{} package is available from the author's website:

@@ -80,7 +80,9 @@ 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.

+are often useful surrogates for batch.  Samples that were processed at

+similar times (e.g., in the same month) can be grouped together in the

+same batch.

 \section{Quality control}

@@ -98,7 +100,7 @@ quality.  The following code chunk makes a histogram of the SNR values

 for the HapMap samples.

 <<snr,fig=TRUE,include=FALSE,width=6, height=4>>=

-open(cnSet$SNR)

+invisible(open(cnSet$SNR))

 snr <- cnSet$SNR[]

 close(cnSet$SNR)

 print(histogram(~snr, panel=function(...){

@@ -120,13 +122,13 @@ print(histogram(~snr, panel=function(...){

 \section{Copy number estimation}

-As described in \cite{Scharpf2010}, the CRLMM-CopyNumber algorithm

+As described in \cite{Scharpf2011}, the CRLMM-CopyNumber algorithm

 fits a linear model to the normalized intensities stratified by the

 diallic genotype call.  The intercept and slope from the linear model

 are both SNP- and batch-specific.  The implementation in the \crlmm{}

 package is encapsulated by the function \Rfunction{crlmmCopynumber}

 that, using the default settings, can be called by passing a single

-object of class \code{CNSet}.  See the appropriate

+object of class \Rclass{CNSet}.  See the appropriate

 preprocessing/genotyping vignette for the construction of an object of

 class \Rclass{CNSet}.

 <<LDS_copynumber,cache=TRUE>>=

@@ -147,7 +149,7 @@ computed for subsets of the markers to reduce the required RAM. Note

 that the value returned by the \Rfunction{crlmmCopynumber} function in

 the above example is \verb+TRUE+.  The reason the function returns

 \verb+TRUE+ in the above example is that the elements of the

-\verb+batchStatistics+ slot have the class \Rclass{ff_matrix}.  Rather

+\verb+batchStatistics+ slot have the class \Rclass{ff\_matrix}.  Rather

 than keep the statistical summaries in memory, the summaries are

 written to files on disk using protocols described in the

 \Rpackage{ff} package. Hence, while the \Robject{cnSet} object itself

@@ -169,9 +171,9 @@ slot of the class \Rclass{CNSet}.  Using the default settings for the

 Note that depends on whether the elements of the \verb+batchStatistics+

 slot are \Robject{ff} objects or ordinary matrices.  In this example,

 the elements of \verb+batchStatistics+ have the class

-\Rclass{ff_matrix}.

+\Rclass{ff\_matrix}.

-<<LDS_copynumber,cache=TRUE>>=

+<<classes>>=

 nms <- ls(batchStatistics(cnSet))

 cls <- rep(NA, length(nms))

 for(i in seq_along(nms)) cls[i] <- class(batchStatistics(cnSet)[[nms[i]]])[1]

@@ -284,6 +286,7 @@ In the following code chunk, we extract estimates of the total copy

 number at nonpolymorphic markers on chromosome X.

 <<npx>>=

+set.seed(123)

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

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

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

@@ -341,9 +344,9 @@ print(bwplot(cn~id,df, panel=function(x,y,...){

 \begin{figure}[t!]

  \centering

  \includegraphics[width=0.5\textwidth]{copynumber-polymorphicX}

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

-  X. crlmm assumes that the median copy number across samples at a

-  given marker on X is 1 for men and 2 for women.}

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

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

 \subsection{A container for raw copy number}

@@ -355,29 +358,31 @@ confidence scores, and the total copy number at each marker is the

 \Rfunction{as} (as illustrated below). Users should note that if the

 \verb+assayData+ elements in the \Rclass{CNSet} instance are

 \Rpackage{ff} objects, the \verb+assayData+ elements of the

-instantiated \Rclass{oligoSnpSet} will also be \Rpackage{ff} object (a

-new \verb+total_cn*.ff+ file will be created in the \verb+ldPath()+

-directory). The advantage of this approach is that the raw copy number

-estimates can be computed for all markers and all samples

-simultaneously without excessive RAM.  The disadvantage is that such a

-step creates additional I/O for reading and writing the raw copy

-number estimates from disk.  For very large data sets (1000+ samples),

-the time required for I/O may outweigh the benefits.  The alternative

-is to perform downstream processing (such as smoothing) on the total

-copy number for subsets of markers and/or samples.  As indicated

-previously, the subset `[' method applied to an object of class

-\Robject{CNSet} will return an object of the same class with simple

-matrices instead of \Rpackage{ff} objects.  Therefore, one can

-instantiate a \Rpackage{oligoSnpSet} object from a \Robject{CNSet}

-object without the I/O overhead for storing the total copy number

-estimates by first subsetting the \Rclass{CNSet} object.  One could

-then smooth the raw copy number estimates and store a summary copy

-number for the interval, thereby substantially reducing the

-dimensionality of the copy number estimates.  Examples of smoothing

-applications such as circular binary segmentation implemented in the

-\Rpackage{DNAcopy} package or a hidden Markov model implemented in the

-\Rpackage{VanillaICE} package are provided in the \verb+SmoothingCN+

-vignette.

+instantiated \Rclass{oligoSnpSet} will also be \Rpackage{ff}-dervied

+objects (a new \verb+total_cn*.ff+ file will be created in the

+\verb+ldPath()+ directory).

+

+%The advantage of this approach is that the raw copy number estimates

+%can be computed for all markers and all samples simultaneously without

+%excessive RAM.  The disadvantage is that such a step creates

+%additional I/O for reading and writing the raw copy number estimates

+%from disk.  For very large data sets (1000+ samples), the time

+%required for I/O may outweigh the benefits.  The alternative is to

+%perform downstream processing (such as smoothing) on the total copy

+%number for subsets of markers and/or samples.  As indicated

+%previously, the subset `[' method applied to an object of class

+%\Robject{CNSet} will return an object of the same class with simple

+%matrices instead of \Rpackage{ff} objects.  Therefore, one can

+%instantiate a \Rpackage{oligoSnpSet} object from a \Robject{CNSet}

+%object without the I/O overhead for storing the total copy number

+%estimates by first subsetting the \Rclass{CNSet} object.  One could

+%then smooth the raw copy number estimates and store a summary copy

+%number for the interval, thereby substantially reducing the

+%dimensionality of the copy number estimates.  Examples of smoothing

+%applications such as circular binary segmentation implemented in the

+%\Rpackage{DNAcopy} package or a hidden Markov model implemented in the

+%\Rpackage{VanillaICE} package are provided in the \verb+SmoothingCN+

+%vignette.

 <<oligoSnpSet>>=

 open(cnSet3)

@@ -388,7 +393,7 @@ class(copyNumber(oligoSet))

 \noindent Note that the raw copy number estimates stored in the

 \Robject{oligoSnpSet} object can be retrieved by the

-\Rfunction{copyNumber} accessor and is equivalanet to that returned by

+\Rfunction{copyNumber} accessor and is equivalent to that returned by

 the \Rfunction{totalCopynumber} function defined over the same row and

 column indices.

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

 \usepackage{graphicx}

 \usepackage{natbib}

 \usepackage{url}

+\usepackage{amsmath,amssymb}

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

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

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

@@ -18,7 +19,7 @@

 \usepackage[margin=1in]{geometry}

 \begin{document}

-\title{Copy number estimation and genotype calling with \Rpackage{crlmm}}

+\title{Copy number estimation with \Rpackage{crlmm}}

 \date{\today}

 \author{Rob Scharpf}

 \maketitle

@@ -27,178 +28,202 @@

 options(continue=" ", width=70)

 @

-%\section{Estimating copy number}

-

-%At present, software for copy number estimation is provided only for the

-%Affymetrix 6.0 platform.

-

 \begin{abstract}

-  This vignette estimates copy number for HapMap samples on the

-  Affymetrix 6.0 platform.  See \citep{Scharpf2010} for details

-  regarding the methodology implemented in \crlmm{}.  In addition, a

-  compendium describing copy number analysis using the \crlmm{}

-  package is available from the author's website:

+

+  Copy number routines in the \crlmm{} package are available for

+  Affymetrix 5.0 and 6.0 platforms, as well as several Illumina

+  platforms.  This vignette assumes that the arrays have already been

+  successfully preprocessed and genotyped as per the instructions in

+  the \verb+AffymetrixPreprocessCN+ and \verb+IlluminaPreprocessCN+

+  vignettes for the Affymetrix and Illumina platforms,

+  respectively. While this vignette uses Affymetrix 6.0 arrays for

+  illustration, the steps at this point are identical for both

+  platforms.  See \citep{Scharpf2010} for details regarding the

+  methodology implemented in \crlmm{} for copy number analysis.  In

+  addition, a compendium describing copy number analysis using the

+  \crlmm{} package is available from the author's website:

   \url{http://www.biostat.jhsph.edu/~rscharpf/crlmmCompendium/index.html}.

 \end{abstract}

-\section{Copy number estimation}

-

-\subsection{Set up}

+\section{Set up}

-<<cdfname>>=

+<<libraries,results=hide>>=

+library(ff)

 library(crlmm)

+library(lattice)

 library(cacheSweave)

-@

-

-Several genotyping platforms are currently supported.  Supported

-platforms must have a corresponding annotation package (installed

-separately) that are listed below.

-

-<<annotationPackages>>=

-pkgs <- annotationPackages()

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

-pkgs

-@

-

-Note that 'pd.genomewidesnp6' and 'genomewidesnp6Crlmm' are both

-annotation packages for the Affymetrix 6.0 platform available on

-Bioconductor, but the 'genomewidesnp6Crlmm' annotation package must be

-used for copy number estimation.  The annotation package is specified

-through the 'cdfName' -- the identifier without the 'Crlmm' postfix.

-In the following code, we specify the cdf name for Affymetrix 6.0,

-provide the complete path to the CEL files, and indicate where

-intermediate files from the copy number estimation are to be saved.

-

-<<setup>>=

-cdfName <- "genomewidesnp6"

-pathToCels <- "/thumper/ctsa/snpmicroarray/hapmap/raw/affy/1m"

 if(getRversion() < "2.13.0"){

 	rpath <- getRversion()

 } else rpath <- "trunk"

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

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

+@

+

+<<ram>>=

+ldPath(outdir)

 setCacheDir(outdir)

+ocProbesets(150e3)

+ocSamples(200)

 @

-All long computations are saved in the output directory

-\Robject{outdir}.  Users should change these variables as appropriate.

-The following code chunk should fail unless the above arguments have

-been set appropriately.

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

+\verb+AffymetrixPreprocessCN+ vignette.

-<<checkSetup>>=

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

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

+<<loadcnset>>=

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

 @

-Processed data from codechunks requiring long computations are saved

-to disk by wrapping function calls in the \Robject{checkExists}

-function.  After running this vignette as a batch job, subsequent