% \VignetteIndexEntry{VanillaICE Vignette}
% \VignetteKeywords{copy number, genotype, SNP}
% \VignettePackage{VanillaICE}
\documentclass{article}
\usepackage{amsmath}
\usepackage{graphicx}

\newcommand{\scscst}{\scriptscriptstyle}
\newcommand{\scst}{\scriptstyle}
\newcommand{\Rpackage}[1]{\textit{#1}}
\newcommand{\Rfunction}[1]{\texttt{#1}}
\newcommand{\Robject}[1]{\texttt{#1}}
\newcommand{\R}{\textsf{R}}
\newcommand{\hmmoptions}{\Robject{HmmOptions}}
\newcommand{\hmmparam}{\Robject{HmmParameter}}

\newcommand{\cne}{\widehat{\text{CN}}}
\newcommand{\gte}{\widehat{\text{GT}}}
\newcommand{\gtehom}{\widehat{\text{HOM}}}
\newcommand{\gtehet}{\widehat{\text{HET}}}
\newcommand{\pgte}{\text{S}_{\widehat{\text{\tiny GT}}}}
\newcommand{\pcne}{\text{S}_{\widehat{\text{\tiny CN}}}}
\newcommand{\pgtehom}{\text{S}_{\widehat{\text{\tiny HOM}}}}
\newcommand{\pgtehet}{\text{S}_{\widehat{\text{\tiny HET}}}}
\newcommand{\thom}{\text{HOM}}
\newcommand{\thet}{\text{HET}}
\newcommand{\bDelta}{\mbox{\boldmath $\Delta$}}
\newcommand{\real}{\mbox{$\mathbb R$}}      % real numbers
\newcommand{\bnu}{\mbox{\boldmath $\nu$}}
\newcommand{\ice}{\Rpackage{VanillaICE}}

\textwidth=6.2in
\textheight=8.5in
\oddsidemargin=.1in
\evensidemargin=.1in
\headheight=-.3in

\begin{document}
\title{\ice{}: Hidden Markov Models for the Assessment of Chromosomal
  Alterations using High-throughput SNP Arrays}
\author{Robert Scharpf}
\maketitle


<<setup, echo=FALSE, results=hide>>=
options(width=69)
@ 

\section{Introduction}

Chromosomal DNA is characterized by variation between individuals at
the level of entire chromosomes (e.g. aneuploidy in which the
chromosome copy number is altered), segmental changes (including
insertions, deletions, inversions, and translocations), and changes to
small genomic regions (including single nucleotide polymorphisms). A
variety of alterations that occur in chromosomal DNA, many of which
can be detected using high density single nucleotide polymorphism
(SNP) microarrays, are linked to normal variation as well as disease
and therefore of particular interest. These include changes in copy
number (deletions and duplications) and genotype (e.g. the occurrence
of regions of homozygosity).  Hidden Markov models (HMM) are
particularly useful for detecting such abnormalities, modeling the
spatial dependence between neighboring SNPs.  Here, we extend previous
approaches that utilize HMM frameworks for inference in high
throughput SNP arrays by integrating copy number, genotype calls, and
the corresponding measures of uncertainty when available.  Using
simulated and real data, we demonstrate how confidence scores control
smoothing in a probabilistic framework.  The goal of this vignette is
to provide a simple interface for fitting HMMs and plotting functions
to help visualize the predicted states alongside the experimental
data.

\section{Simple Usage}
\label{sec:simpleUsage}

\subsection{Vanilla HMM }

<<packages, results=hide>>=
library(VanillaICE)
@ 

Before fitting the HMM, the data must be organized into one of the
following classes for high-throughput SNP data:

\begin{itemize}
\item \Robject{SnpCallSet} (genotype calls)
\item \Robject{SnpCopyNumberSet} (copy number estimates)
\item \Robject{oligoSnpSet} (genotype calls and copy number estimates)
\item \Robject{RatioSnpSet} (genotype calls and ratios of A and B allele intensities)
\end{itemize}

When pre-processing Affymetrix SNP chips with the \R{} package
\Rpackage{oligo}, an object of one of the above classes is created and
can be used directly with the HMMs described in this vignette.  These
classes can also be created from Illumina data as described in the
IlluminaHowTo vignette.  The simulated data provided with this package
is an instance of class \Robject{oligoSnpSet}

<<data>>=
data(chromosome1)
annotation(chromosome1)
chromosome1
@ 

Before deciding whether to fit a HMM, plot the data:

<<plotSnp,fig=TRUE,width=8, height=5>>=
gp <- plotSnp(chromosome1)
show(gp)
@ 

The HMM assumes that the copy number estimates, conditional on the
hidden state, are approximately Gaussian.  A log transformation will
be performed automatically provided that the copy number estimates are
all positive.  See the documentation pages in the R package \ice{} for
more information about the \Robject{chromosome1} example dataset.
Fitting a hidden Markov model requires the following components:

\begin{itemize}
  \item the hidden states
  \item the emission probabilities
  \item transition probabilities
\end{itemize}

Here, we illustrate how one may specify the hidden states and
calculate the emission and transition probabilities below.  We first
create an instance of the class \hmmoptions{} that contains the hidden
states and, if copy number estimates are available, the location
parameters for the Gaussian distribution.  Though we will model the
logarithm of the copy number estimates as approximately normally
distributed, it is more natural to specify the location parameters on
the original scale (e.g., 1 copy is deletion (D), 2 copies is normal
(N), etc.).  Emission probabilities for the genotype calls are a
homozygous genotype call (\Robject{probHomCall} argument).  Again, the
ordering must correspond to the order of the hidden states.  computed
from the supplied probabilities of Note that the ordering of the
numeric vector for copyNumber.location must correspond to the ordering
of the states.  A number of validity checks are performed on this
object and it is a good idea to check whether the instantiated object
is valid.

<<options>>=
options <- new("HmmOptions",
	       states=c("D", "N", "L", "A"),
	       snpset=chromosome1,
	       copyNumber.location=c(1, 2, 2, 3),
	       probHomCall=c(0.99, 0.7, 0.99, 0.7))
validObject(options)
@ 

From an object of class \hmmoptions{}, we can compute the emission
probabilities.  The emission probabilities are returned on the log
scale.  Conditional on the hidden state, we assume that the copy
number and genotype are independent.  Therefore, the emission
probabilities for an HMM that models the copy number and genotypes
jointly are computed by adding the emission probabilities for copy
number and genotype.

<<HmmParameters>>=
params <- new("HmmParameter", 
	      states=states(options),
	      initialStateProbability=0.99)
cn.emission <- copyNumber.emission(options)
gt.emission <- calls.emission(options)
emission(params) <- cn.emission + gt.emission ##log scale
@ 

<<genomicDistance>>=
genomicDistance(params) <- exp(-2 * physicalDistance(options)/(100*1e6))
transitionScale(params) <- scaleTransitionProbability(options)
class(params)
@ 

We may then fit the HMM by

<<fit, eval=FALSE>>=
fit <- hmm(options, params)
@ 

<<load,echo=FALSE>>=
data(fit)
@ 

<<>>=
fit
breakpoints(fit)
@ 

Summary statistics for the breakpoints (more useful when multiple
chromosomes and samples are in the \Robject{fit} object)

<<summary>>=
summary(fit)
@ 


See \cite{Scharpf2007a} for a more complete description of the
simulated dataset and the features detected by this HMM.  We may plot
the data along with the predictions as follows:

<<vanillaPlot,fig=TRUE,width=8,height=5>>=
gp <- plotSnp(snpset(options), fit)
gp$abline.v <- TRUE ##plots vertical dashed lines
show(gp)
@ 

\subsection{Integrating Confidence Estimates (ICE)}

In this section, we illustrate how one may fit an HMM that
incorporates confidence esimates of the SNP-level summaries for
genotype calls and copy number.  Confidence scores (inverse of
standard errors) are available for this object (see Section
\ref{sec:confidenceScores} for how confidence scores were derived).

This information is incorporated into the HMM emission
probabilities. Probably the easiest way to do that is recreate the
options object, and then recalculate the emission probabilities.


<<params>>=
options <- new("HmmOptions", 
	       snpset=chromosome1, 
	       states=c("D", "N", "L", "A"), 
	       copyNumber.location=c(1, 2, 2, 3),
	       copyNumber.ICE=TRUE,
	       probHomCall=c(0.99, 0.75, 0.99, 0.75))
cn.emission <- copyNumber.emission(options)
emission(params) <- cn.emission + gt.emission
fit.ice <- hmm(options, params)
@ 

We may also incorporate the confidence scores for the genotype calls
by specifying \Robject{calls.ICE=TRUE}. This feature is only available
for the Affymetrix 100k and 500k platforms. The slot
\Robject{probHomCall} stores user-specified probabilities of P(call is
AA or BB | state is LOH) and P(call is AA or BB | state is normal).
These probabilities must be specified in this order.  The emission
probabilities for the genotype calls will only be calculated for the
states LOH (LOH is defined as a stretch of homozygous genotype calls
longer than what one would expect by chance) and Normal refers to
typical ratios of heterozygous to homozygous genotype calls. The slot
\Robject{term5} contains user-specified probabilities for the P(true
genotype is HET | genotype call is AB, hidden state is Normal) and
P(true genotype is HET | genotype call is AA or BB, hidden state is
Normal), respectively.  Default values are provided when not
specified, as the following example illustrates.

<<params2, eval=FALSE>>=
options <- new("HmmOptions", 
	       snpset=chromosome1, 
	       states=c("D", "N", "L", "A"), 
	       copyNumber.location=c(1, 2, 2, 3),
	       copyNumber.ICE=TRUE, 
	       calls.ICE=TRUE,
	       probHomCall=c(0.99, 0.75))
params <- new("HmmParameter", 
	      states=states(options),
	      initialStateProbability=0.99)
cn.emission <- copyNumber.emission(options)
genomicDistance(params) <- exp(-2 * physicalDistance(options)/(100*1e6))
transitionScale(params) <- scaleTransitionProbability(options)
gt.emit <- calls.emission(options)
gt.emission <- array(NA, dim(cn.emission))
gt.emission[, , 1:2] <- gt.emit
gt.emission[, , 3:4] <- gt.emit
emission(params) <- cn.emission + gt.emission
fit.ice <- hmm(options, params)
@ 


<<icePlot, fig=TRUE, width=8, height=5>>=
gp <- plotSnp(snpset(options), fit.ice)
gp$abline.v <- TRUE
gp$col.predict <- c("darkblue", "white", "yellow", "orange")
show(gp)
legend(-0.05, 10, fill=gp$col.predict[c(1, 3, 4)],
       legend=c("< 2 copies", "copy-neutral LOH", "> 2 copies"),
       bty="n", title="predicted states")
legend(0, 0.8, legend="predictions", bty="n", cex=0.8, adj=0)
@

Note that the ICE HMM correctly identifies the simulated normal
segments in features B and C (the normal segments were simulated to
have high confidence scores). Additionally, the ICE HMM detects the
micro-amplification in region E (also simulated to have high
confidence scores).

\subsection{Confidence scores}\label{sec:confidenceScores}

\paragraph{Confidence scores for genotype calls}

We suggest using the \Rfunction{CRLMM} algorithm \cite{Carvalho2007}
for genotype calls.  \Rfunction{CRLMM} (in the R package
\Rpackage{oligo}) provides confidence scores ($\pgte$) of the genotype
estimates ($\gte$).  From 269 HapMap samples assayed on the Affymetrix
50k Xba and Hind chips, we have a gold standard of the true genotype
defined by the consensus of the HapMap centers.  We use kernal based
density estimates to obtain

{\scriptsize
\begin{eqnarray}
\label{ingo2}
f\left\{\ \pgtehom\ |\ \gtehom,\thom\ \right\},\
f\left\{\ \pgtehom\ |\ \gtehom,\thet\ \right\},\
f\left\{\ \pgtehet\ |\ \gtehet,\thom\ \right\},\ \mbox{~ and~}
f\left\{\ \pgtehet\ |\ \gtehet,\thet\ \right\}
\end{eqnarray}
}

\noindent separately for the Xba and Hind 50k chips. The first term in
(\ref{ingo2}), for example, denotes the density of the scores when the
genotype is correctly called homozygous ($\gtehom$) and the true
genotype is homozygous ($\thom$). See \cite{Scharpf2007a} for a more
complete description of the methods. The data needed to estimate these
densities is stored in the experiment data package
\Rpackage{callsConfidence}.  \Rpackage{callsConfidence} is available
from the author's website.
%\texttt{http://biostat.jhsph.edu/~rscharpf/software/index.html}.

\paragraph{Confidence scores for copy number estimates}

To illustrate how standard errors of the copy number estimate could be
integrated in the HMM, the R object \Robject{chromosome1} contains
standard errors simulated from a shifted Gamma: $\mbox{Gamma}(1, 2) +
0.3$, where 1 is the shape parameter and 2 is the rate parameter.  To
ascertain the effect of qualitatively high confidence scores on the
ICE HMM, we scaled a robust estimate of the copy number standard
deviation by $\frac{1}{2}$. Similarly, to simulate less precise $\cne$
we scaled $\epsilon$ by 2.  For more detailed information about how
the data in the \Robject{chromosome1} was generated, see the
documentation for this object in the R package \ice.

\section{The HmmParameter class}

  
An instance of the class is created by the method \Rfunction{new}:
  
<<hmmParameters, eval=FALSE>>=
new("HmmParameter")
@

The object \Robject{params} contains all of the parameters needed for
fitting the HMM, including an estimate of the genomic distance between
SNPs (used for calculating SNP-specific transition probabilities),
emission probabilities (slot: emission), and initial state
probabilities.  


\paragraph{Emission probabilities.}
The emission probabilities are stored as an array in the
\Robject{params} object.  The emission probability array has dimension
R $\times$ C $\times$ S, where $S$ is the number of hidden states, $R$
is the number of rows (SNPs), and $C$ is the number of samples.  One
may use $[$ to subset object of class \Robject{HmmParameter}.

<<>>=
params[5, 1, ]
@ 


\paragraph{Transition probabilities.}
The probability of remaining in the same state, $P(S_t = S_{t+1})$
(the diagonal of the transition probability matrix) is a function of
the distance (d) between SNPs: $e^{-2d(100 * 1e6)}$.  This value is
stored in the slot \Robject{tau} of the \Robject{params} object.  The
probability of leaving a state is $\epsilon$, where $\epsilon = 1 -
P(S_t = S_{t+1})$.  The $\epsilon$ is split among $S-1$ states.  By
default, the probability of transitioning from an altered state back
to the normal state is twice as likely as the probability of
transitioning between two altered states.  The weights for $\epsilon$
are provided in the \texttt{tau.scale} matrix in the R object
\Robject{params}

<<tau.scale>>=
transitionScale(params)
@

\noindent and can be adjusted by the \texttt{SCALE} argument.  For
illustration, one could make the probability of transitioning from an
altered state to a normal state 10 times as likely as the probability
of transitioning between two altered states by the following commmand:

<<>>=
transitionScale(params) <- scaleTransitionProbability(params, SCALE=10)
transitionScale(params)
@ 

\section{The HmmOptions Class}

To be completed \ldots

\section{The HmmPredict Class}

The output from the HMM is an instance of the \Robject{HmmPredict}
class and contains the predicted states as well as the breakpoints for
the different states.

<<hmmOutput>>=
fit
@ 

The breakpoints are provided as a \Robject{data.frame}:

<<breaks>>=
breaks <- breakpoints(fit)
breaks <- breaks[breaks[, "state"] != "N", ]
@ 

One may order the altered states from biggest to smallest for each
chromosome as follows:

<<orderBreaks>>=
breaks[order(breaks[, "chr"], breaks[, "size"], decreasing=TRUE), ]
@ 

The \Robject{summary} method returns a list where each element in the
list provides statistics for an altered states.  For each chromosome,
the mean, median, and standard deviation of the size of the features
and number of SNPs involved are reported. For instance, in this HMM
the altered states were loss of heterozygosity (L), amplification of
copy number (A), and deletion of copy number (D). Because we only have
one sample and one chromosome in the \Robject{fit} example, the
\Robject{summary} method is not that useful:

<<summary>>=
summary(fit)
@ 

\section{HMMs for different classes of data}

\subsection{Copy number}

The method \Robject{hmm} has a different set of underlying hidden
states depending on whether copy number estimates, genotype calls, or
both are available.  When only copy number estimates are available,
the hidden states (for autosomes) are hemizygous or homozygous
deletion (one or fewer copies), normal (two copies), and amplification
(three or more copies).  The corresponding data class is
{SnpCopyNumberSet}.  To illustrate, we convert the
\Robject{chromosome1} example to an object of this class and fit the
HMM.

<<copyNumberHmm>>=
chr1.cn <- as(chromosome1, "SnpCopyNumberSet")
options <- new("HmmOptions",
	       snpset=chr1.cn,
	       states=c("D", "N", "A"),
	       copyNumber.location=1:3)
params.cn <- new("HmmParameter",  
		 states=c("D", "N", "A"))
emission(params.cn) <- copyNumber.emission(options)
transitionScale(params.cn) <- scaleTransitionProbability(options)
genomicDistance(params.cn) <- exp(-2*physicalDistance(options)/(100*1e6))
fit.cn <- hmm(options, params.cn)
breakpoints(fit.cn)
graph.par <- plotSnp(snpset(options), fit.cn)
graph.par$abline.v <- FALSE
@ 

<<copyNumberFigure, fig=TRUE, include=TRUE>>=
print(graph.par)
legend(0, 0.8, fill=graph.par$col.predict[c(1, 3)], 
       legend=c("< 2 copes", "> 2 copies"),
       bty="n", cex=0.8)
@ 

<<echo=FALSE>>=
rm(chr1.cn, fit.cn, params.cn)
@ 

\subsection{Genotype calls}

When only genotype calls are available, the hidden states are loss and
retention (ret) of heterozygosity.  We define \textit{loss} to be a
sequence of homozygous SNPs longer than what we would expect to
observe by chance. Note that many long stretches of homozygosity may
occur as a result of a population sharing a common underlying
haplotype structure; loss predictions from an HMM fit to an indvidual
do not necessarily reflect the 'loss' of an allele in that individual.
For illustration, we convert the \Robject{chromosome1} example to an
object of class \Robject{HmmSnpCallSet} and refit the HMM.

<<chr1Calls>>=
chr1.calls <- as(chromosome1, "SnpCallSet")
options <- new("HmmOptions",
	       snpset=chr1.calls,
	       states=c("L", "N"),
	       probHomCall=c(0.99, 0.7))
params.calls <- new("HmmParameter", 
		    states=states(options))
transitionScale(params.calls) <- scaleTransitionProbability(options)
genomicDistance(params.calls) <- exp(-2*physicalDistance(options)/(100*1e6))
emission(params.calls) <- calls.emission(options)

fit.calls <- hmm(options, params.calls)
breakpoints(fit.calls)
gp <- plotSnp(snpset(options), fit.calls)
gp$col.predict <- c("black", "white")
gp$ylim <- c(-0.5, 1)
gp$add.centromere <- FALSE
gp$abline.v <- TRUE
gp$cytoband.ycoords <- c(-0.45, -0.4)
gp$hmm.ycoords <- c(-0.2, -0.05)
@ 

<<callsFigure, fig=TRUE, include=TRUE>>=
show(gp)
legend(0, -0.1, legend="LOH", fill="black", 
       title="predictions", bty="o", cex=0.8)
@ 

\subsection{Genotype calls and copy number}

Section \ref{sec:simpleUsage} illustrates how one may fit the HMM to
objects of class \Robject{oligoSnpSet}.

More documentation about the classes can be found in the documentation
for the R package \ice. 

\section{Session Information}

The version number of R and packages loaded for generating the
vignette were:

<<echo=FALSE, results=tex>>=
toLatex(sessionInfo())
@ 

\bibliography{ice}{}
\bibliographystyle{plain}

\end{document}
