deleted file mode 100755

@@ -1,287 +0,0 @@

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

-%\VignetteDepends{crlmm}

-%\VignetteKeywords{crlmm, illumina}

-%\VignettePackage{crlmm}

-\documentclass{article}

-\usepackage{graphicx}

-\usepackage{natbib}

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

-\newcommand{\crlmm}{\Rpackage{crlmm}}

-\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{\ff}{\Rpackage{ff}}

-

-\begin{document}

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

-  with Illumina platforms}

-

-\date{\today}

-

-\author{Rob Scharpf}

-\maketitle

-

-

-\begin{abstract}

-  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{Setup}

-

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

-library(crlmm)

-options(width=70)

-options(continue=" ")

-@

-

-%\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.

-The following codechunk declares a directory for saving \Robject{ff}

-files that will contain the normalized intensities and the genotype

-calls.

-

-<<ldpath,results=hide>>=

-library(ff)

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

-@

-

-We will also store cached computations in the directory \verb+outdir+.

-

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

-library(cacheSweave)

-setCacheDir(outdir)

-@

-

-We declare that \crlmm{} should process 150,000 markers at a time

-(when possible) and/or 500 samples at a time.  As our example dataset

-in this vignette contains fewer than 500 samples, all samples will be

-processed simultaneously.

-

-<<ram>>=

-ocProbesets(150e3)

-ocSamples(500)

-@

-

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

-

-<<datadir>>=

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

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

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

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

-@

-

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

-

-<<cdfname>>=

-cdfName <- "human370v1c"

-@

-

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

-

-<<batch>>=

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

-@

-

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

-

-<<cnset>>=

-print(cnSet)

-@

-

-Note that the above object does not yet contain any processed data

-(only \verb+NA+'s) and several \verb+.ff+ files now appear in the

-\verb+outdir+.  Again, these files should not be removed.

-

-<<listff>>=

-list.files(outdir, pattern=".ff")[1:5]

-@

-

-Finally, note that the elements of the \verb+assayData+ slot are

-\Robject{ff} objects and not matrices. For the most part, the

-\emph{appearance} that the data is stored in memory is preserved.

-

-

-<<assaydataelements>>=

-sapply(assayData(cnSet), function(x) class(x)[1])

-@

-

-\section{Preprocessing}

-

-The raw intensities from the Infinium IDAT files are read and

-normalized using the function \Rfunction{preprocessInf}.  The function

-\Rfunction{preprocessInf} returns a \Robject{ff} object containing the

-parameters for the mixture model used by the CRLMM genotyping

-algorithm.

-

-<<preprocess,cache=TRUE>>=

-mixtureParams <- preprocessInf(cnSet=cnSet, sampleSheet=samplesheet, arrayNames=arrayNames, arrayInfoColNames=arrayInfo)

-invisible(open(mixtureParams))

-str(mixtureParams[])

-invisible(close(mixtureParams))

-@

-

-Note that the normalized intensities for the A and B alleles are no

-longer \verb+NA+s and can now be accessed and inspected using the

-methods \Rfunction{A} and \Rfunction{B}, respectively.

-

-<<intensities>>=

-invisible(open(A(cnSet)))

-invisible(open(B(cnSet)))

-as.matrix(A(cnSet)[1:5, 1:5])

-as.matrix(B(cnSet)[1:5, 1:5])

-invisible(close(A(cnSet)))

-invisible(close(B(cnSet)))

-@

-

-

-\section{Genotyping}

-

-

-CRLMM genotype calls and confidence scores are estimated using the

-function \Rfunction{genotypeInf}.

-

-<<genotype,cache=TRUE>>=

-updated <- genotypeInf(cnSet, mixtureParams=mixtureParams)

-print(updated)

-@

-

-The posterior probabilities for the

-genotype calls in the \verb+callProbability+ element of the

-\verb+assayData+ are stored as integers to reduce the file size on

-disk.  However, the scores can be easily transformed back to the

-probability scale using the \Rfunction{i2p} function as illustrated in

-the following code chunk.

-

-<<callprobs>>=

-invisible(open(snpCallProbability(cnSet)))

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

-i2p(callProbs)

-invisible(close(snpCallProbability(cnSet)))

-@

-

-%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.

-

-\textbf{Wrapper:} The construction of a \Rclass{CNSet} object,

-preprocessing, and genotype calling are wrapped into a convenience

-function called \Rfunction{genotype.Illumina}.

-

-<<genotype.Illumina,cache=TRUE,results=hide>>=

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

-			    arrayNames=arrayNames,

-			    arrayInfoColNames=arrayInfo,

-			    cdfName="human370v1c",

-			    batch=batch)

-@

-

-<<checkIdenticalCalls>>=

-invisible(open(calls(cnSet)))

-invisible(open(calls(cnSet2)))

-snp.index <- which(isSnp(cnSet))

-identical(calls(cnSet)[snp.index, 1:20], calls(cnSet2)[snp.index, 1:20])

-invisible(close(calls(cnSet)))

-invisible(close(calls(cnSet2)))

-@

-

-To fully remove the data associated with the \Robject{cnSet2} object,

-one should use the \Rfunction{delete} function in the \Rpackage{ff}

-package followed by the \Rfunction{rm} function.  The following code

-is not evaluated is it would change the results of the cached

-computations in the previous code chunk.

-

-<<delete, eval=FALSE>>=

-lapply(assayData(cnSet2), delete)

-lapply(batchStatistics(cnSet2), delete)

-delete(cnSet2$gender)

-delete(cnSet2$SNR)

-delete(cnSet2$SKW)

-rm(cnSet2)

-@

-

-\end{document}

@@ -36,7 +36,7 @@

   to the copy number vignette.

 \end{abstract}

-\section{Infrastructure}

+\section{Setup}

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

 library(crlmm)

@@ -44,29 +44,16 @@ options(width=70)

 options(continue=" ")

 @

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

-

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

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

+%\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.

+The following codechunk declares a directory for saving \Robject{ff}

+files that will contain the normalized intensities and the genotype

+calls.

 <<ldpath,results=hide>>=

 library(ff)

@@ -78,41 +65,23 @@ ldPath(outdir)

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

 @

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

-ocSamples()

-ocProbesets()

-@

+We will also store cached computations in the directory \verb+outdir+.

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

 library(cacheSweave)

 setCacheDir(outdir)

 @

+We declare that \crlmm{} should process 150,000 markers at a time

+(when possible) and/or 500 samples at a time.  As our example dataset

+in this vignette contains fewer than 500 samples, all samples will be

+processed simultaneously.

+

+<<ram>>=

+ocProbesets(150e3)

+ocSamples(500)

+@

+

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

@@ -134,7 +103,7 @@ 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>>=

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

 @

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

-

-Alternatively, total copy number can be obtained by

-<<totalCopynumber.snps>>=

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

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

+<<cnset>>=

+print(cnSet)

 @

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

+Note that the above object does not yet contain any processed data

+(only \verb+NA+'s) and several \verb+.ff+ files now appear in the

+\verb+outdir+.  Again, these files should not be removed.

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

+<<listff>>=

+list.files(outdir, pattern=".ff")[1:5]

 @

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

-@

-

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

-

-<<copynumberObject>>=

-marker.index <- which(chromosome(cnSet) <= 22)## & isSnp(cnSet))

-sample.index <- 1:5 ## first five samples

-invisible(open(cnSet))

-copynumberSet <- as(cnSet[marker.index, 1:5], "CopyNumberSet")

-invisible(close(cnSet))

-copynumberSet <- copynumberSet[order(chromosome(copynumberSet), position(copynumberSet)), ]

-indices <- split(1:nrow(copynumberSet), chromosome(copynumberSet))

-dup.index <- unlist(sapply(indices, function(i, position)  i[duplicated(position[i])], position=position(copynumberSet)))

-if(length(dup.index) > 0) copynumberSet <- copynumberSet[-dup.index, ]