%\VignetteIndexEntry{Infrastructure for copy number analysis}
%\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

\begin{abstract}

  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 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).  This object was instantiated by running
  the \verb+AffyGW+ vignette.

\end{abstract}

<<libraries,results=hide>>=
library(ff)
library(crlmm)
@

\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 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(100e3)
ocSamples(200)
@

\section{The \Rclass{CNSet} container}

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

An object of class \Rclass{CNSet} can be instantiated by one of two
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 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
  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{itemize}

There are important differences in the underlying data representation
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}

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>>=
outdir <- paste("/local_data/r00/crlmm/", getRversion(), "/infrastructure", sep="")
ldPath(outdir)
if(!file.exists(outdir)) dir.create(outdir)
@


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+AffyGW+ vignette.

<<instantiateToyExample,cache=TRUE>>=
exampleSet <- genotype(celfiles, batch=rep("1", 3), cdfName="genomewidesnp6")
@

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 exported function
\Rfunction{constructAffy} directly using the \Rfunction{:::} operator.

<<constructAffy,eval=FALSE>>=
tmp <- 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
using accessors to extract data from the class.  For example, the
following 2 methods would extract identical results, with the latter
being much more efficient (extra parentheses are added to the second
operation to emphasize the order of operations):

<<orderOfOperations>>=
##inefficient
##invisible(open(cnSet))
A(exampleSet[1:5, ])
## preferred
(A(exampleSet))[1:5, ]
@


\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>>=
library(Biobase)
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(exampleSet) == 1)
@

\subsubsection{\texttt{assayData}}

The \verb+assayData+ elements are of the class
\Rclass{ff\_matrix}/\Rclass{ffdf} or \Rclass{matrix}, depending on how
the \Rclass{CNSet} object was instantiated.  Elements in the
\verb+assayData+ environment can be listed using the \Rfunction{ls}
function.

<<assayData>>=
ls(assayData(exampleSet))
@

The normalized intensities for the A and B alleles have names
\verb+alleleA+ and \verb+alleleB+ and can be accessed by the methods
\Rfunction{A} and \Rfunction{B}, respectively.  Genotype calls and
confidence scores can be accessed by \Rfunction{snpCall} and
\Rfunction{snpCallProbability}, respectively.  Note that the
confidence scores are represented as integers to reduce the filesize,
but can be tranlated to the probability scale using the \R{} function
\Rfunction{i2p}. For example,

<<i2p>>=
scores <- as.matrix(snpCallProbability(exampleSet)[1:5, 1:2])
i2p(scores)
@

Note that for the Affymetrix 6.0 platform the assay data elements each
have a row dimension corresponding to the total number of polymorphic
and nonpolymorphic markers interrogated by the Affymetrix 6.0
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.

\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
  behaviour across conditions and are unrelated to the biological or
  scientific variables in a study}.  The \verb+batchStatistics+ slot
is an environment used to store SNP- and batch-specific summaries,
such as the sufficient statistics for the genotype clusters and the
linear model parameters used for copy number estimation.  The
\verb+batch+ slot is used to store the 'batch name' for each
array. For small studies in which the samples were processed at
similar times (e.g., within a month), all the samples can be
considered to be in the same batch.  For large studies in which the
samples were processed over several months, users should the scan date
of the array or the chemistry plate are useful surrogates. The only
constraint on the \verb+batch+ variable is that it must be a character
vector that is the same length as the number of samples to be
processed.  The \verb+batch+ is specified as an argument to the \R{}
functions \Rfunction{constructInf} and \Rfunction{genotype} that
instantiate \Rclass{CNSet} objects for the Illumina and Affymetrix
platforms, respectively.  The \Rfunction{batch} function can be used
to access the \verb+batch+ information on the samples as in the
following example.

<<batchAccessor>>=
batch(exampleSet)
@

For the \verb+batchStatistics+ slot, the elements in the environment
have the class \Rclass{ff\_matrix}/\Rclass{ffdf} or \Rclass{matrix},
depending on how the \Rclass{CNSet} object was instantiated.  The
dimension of each element is the number of markers (SNPs +
nonpolymorphic markers) $\times$ the number of batches. The names of
the elements in the environment can be list using the \R{} function
\Rfunction{ls}.

<<batchStatistics>>=
ls(batchStatistics(exampleSet))
@

Currently, the batch-specific summaries are stored to allow some
flexibility in the choice of downstream analyses of copy number and
visual assessments of model fit.  Documentation for such applications
will be expanded in future versions of \crlmm{}, and are currently not
intended to be accessed directly by the user.

%The \Robject{CNSet} class also contains the slot
%\Robject{batchStatistics} that contains batch-specific summaries
%needed for copy number estimation. In particular, each element is a
%matrix (or an ff object) with R rows and C columns, correspoinding to
%R markers and C batches.  The summaries includes the within genotype
%cluster medians and median absolute deviations (mads), but also
%parameters estimated from the linear model.  (For unobserved
%genotypes, the medians are imputed and the variance is obtained the
%median variance (across markers) within a batch. ) The elements of the
%slot can be listed as follows.

\subsubsection{\texttt{phenoData}}

Sample-level summaries obtained during the preprocessing/genotyping
steps include skew (SKW), the signal to noise ratio (SNR), and gender
(1=male, 2=female) are stored in the \verb+phenoData+ slot.  As for
other \Rclass{eSet} extensions, the \verb+$+ method can be used to
extract these summaries.  For \Rclass{CNSet} objects generated by
approach (1), these elements are of the class
\Rclass{ff\_vector}.

<<phenodataAccessors>>=
varLabels(exampleSet)
class(exampleSet$gender)
invisible(open(exampleSet$gender))
exampleSet$gender
@

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

<<vector>>=
c("male", "female")[exampleSet$gender[]]
invisible(close(exampleSet$gender))
@




\subsubsection{\texttt{protocolData}}

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

<<protocolData>>=
varLabels(protocolData(exampleSet))
protocolData(exampleSet)$ScanDate
@

\section{Trouble shooting with a HapMap example}

This section uses an object of class \Rclass{CNSet} instantiated by
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.

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

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).  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 <- 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)
missing.cb <- is.na(cb)
(nmissing.ca <- sum(missing.ca))
(nmissing.cb <- sum(missing.cb))
identical(nmissing.ca, nmissing.cb)
@

If \Robject{nmissing.ca} is nonzero, check the genotype confidence
scores provided by the crlmm genotyping algorithm against the
threshold specified in \Robject{crlmmCopynumber}.

<<NAconfidenceScores>>=
if(nmissing.ca > 0){
	##invisible(open(snpCallProbability(cnSet)))
	gt.conf <- i2p(snpCallProbability(cnSet)[autosome.index, sample.index])
	##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)
	stopifnot(identical(nmissingBecauseOfGtThr, nmissing.ca))
	## or calculate the proportion of missing effected by low crlmm confidence
	length(index.allbelow) * length(sample.index)/nmissing.ca
}
@

One could inspect the cluster plots for the 'low confidence' calls.

<<clusterPlotsLowConfidence>>=
## TODO
@

We repeat the above check for missing values at polymorphic loci on
chromosome X.  In this case, we compare the \Robject{rowSums} of the
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 <- 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)
(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
@

From the above codechunk, we see that \Sexpr{length(index)} SNPs have
NAs for all the samples. Next, we tally the number of NAs for
polymorphic markers on chromosome X that are below the confidence
threshold.  For the HapMap analysis, all of the missing values arose
from SNPs in which either the men or the women had confidence scores
that were all below the threshold.

<<NAcrlmmConfidenceScore>>=
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
and estimates are available at most markers.

<<NAnonpolymorphic.autosome>>=
np.index <- which(!isSnp(cnSet) & chromosome(cnSet)==23)
ca.F <- CA(cnSet, i=np.index, j=F)
ca.M <- CA(cnSet, i=np.index, j=M)
## NAs for one marker
ca.F <- ca.F[-match("CN_974939", rownames(ca.F)), ]
ca.M <- ca.M[-match("CN_974939", rownames(ca.M)), ]
sum(is.na(ca.F))
sum(is.na(ca.M))
@

%TODO: marker CN\_974939 has NAs for the normalized intensities.  This
%is because CN\_974939 is not in the \Robject{npProbesFid} file in
%\Rpackage{genomewidesnp6Crlmm}.  The \Robject{npProbesFid} file should
%be updated in the next \Rpackage{genomewidesnp6Crlmm} release.

In total, there were \Sexpr{length(missing.snp.index)} polymorphic
markers on chromosome X for which copy number estimates are not
available.  Lowering the confidence threshold would permit estimation
of copy number at most of these loci.  A confidence threshold is
included as a parameter for the copy number estimation as an approach
to reduce the sensitivity of genotype-specific summary statistics,
such as the within-genotype median, to intensities from samples that
do not clearly fall into one of the biallelic genotype clusters. There
are drawbacks to this approach, including variance estimates that can
be a bit optimistic at some loci.  More direct approaches for outlier
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{}.

<<close>>=
invisible(close(cnSet))
invisible(close(cnSet$gender))
@


\section{Session information}
<<sessionInfo, results=tex>>=
toLatex(sessionInfo())
@

%\section*{References}
%
%%\begin{bibliography}
%  \bibliographystyle{plain}
%  \bibliography{refs}
%\end{bibliography}


\end{document}