%\VignetteIndexEntry{cape}
%\VignetteDepends{}
%\VignetteKeywords{cape epistasis pleiotropy analysis}
%\VignettePackage{cape}

\documentclass[11pt]{article}

\usepackage{amsmath}
\usepackage[authoryear,round]{natbib}
\usepackage{hyperref}
\usepackage{color}
\usepackage{rotating}
\usepackage[none]{hyphenat}
\usepackage[english]{babel}
\usepackage{blindtext}
\usepackage{array}
\usepackage{ragged2e}
\usepackage{microtype}
\usepackage{caption}

\newcolumntype{P}[1]{>{\RaggedRight\hspace{0pt}}p{#1}}

\setlength{\textheight}{8.5in}
\setlength{\textwidth}{6in}
\setlength{\topmargin}{-0.25in}
\setlength{\oddsidemargin}{0.25in}
\setlength{\evensidemargin}{0.25in}
\setlength\parindent{0pt}
\usepackage{Sweave}


\setkeys{Gin}{width=\textwidth} %% <- change width of figures
% Try to get the R code from running into the margin
%<<echo=FALSE>>=
%options(width=87, digits=3, scipen=4)
%set.seed(61777369)
%@ 

\begin{document}
\setkeys{Gin}{width=0.99\textwidth}
\renewcommand{\arraystretch}{1.5}


\title{\bf cape: A package for the combined analysis of epistasis and pleiotropy}

\author{Anna L. Tyler$^1$, Wei Lu$^{1,2}$, Justin J. Hendrick$^1$, Vivek M. Philip$^1$, and Gregory W. Carter$^1$}

\maketitle

\begin{center}
$^1$The Jackson Laboratory, Bar Harbor, ME, 04609
$^2$Department of Electrical and Computer Engineering, Duke University, Durham, NC 27708


\end{center}

\tableofcontents

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Abstract}
\addcontentsline{toc}{section}{Abstract}

Here we present an R package for the Combined Analysis of Epistasis and 
Pleiotropy, or \texttt{cape}. This package implements a method, originally described
in \cite{Carter:2012fd}, that infers directed interaction networks between genetic 
variants for predicting the influence of genetic perturbations on phenotypes. This 
method takes advantage of complementary information in partially pleiotropic genetic variants 
to resolve directional influences between variants that interact epistatically. \texttt{cape}
can be applied to a variety of genetic variants, such as single nucleotide polymorphisms
(SNPs), copy number variations (CNVs) or structural variations (SVs). Here we 
demonstrate the functionality of \texttt{cape} by inferring a predictive network between
quantitative trait loci (QTL) in a cross between the non-obese, non-diabetic (NON) mouse
and the New Zealand obese (NZO) mouse \citep{Reifsnyder:2000is}.

\section*{Loading Data}
\addcontentsline{toc}{section}{Loading Data}

For the purposes of demonstration, we will reanalyze a data set described in 
\cite{Reifsnyder:2000is}. This data set was established to find quantitative trait 
loci (QTL) for obesity and other risk factors of type II diabetes in a reciprocal 
back-cross of non-obese non-diabetic NON/Lt mice and diabetes-prone, New 
Zealand obese (NZO/HILt) mice. The study found multiple main-effect QTL 
influencing phenotypes associated with diabetes and obesity as well as multiple 
epistatic interactions. In addition, maternal environment (i.e. whether the mother 
was obese) was found to interact with several markers and epistatic pairs 
to influence the risk of obesity and diabetes of the offspring. The complex nature 
of diabetes and obesity, along with their complex and polygenic inheritance 
patterns, make this data set ideal for an analysis of epistasis and pleiotropy.\\

Included in this dataset are 204 male mice genotyped at 85 markers
across the genome. The phenotypes included are the body weight (g), 
insulin levels (ng/mL), and plasma glucose levels (mg/dL), all measured 
at age 24 weeks. In addition, there is a variable called ``mom'' indicating 
whether the mother of each mouse was normal weight (0) or obese (1). 
After installing \texttt{cape}, load the package and the dataset type the 
following in the R command line.

<<>>=
library(cape)
data(obesity.cross)
@

To load your own data, use the function \texttt{{\color{blue}read.population()}}. 
This function is similar to the function \texttt{{\color{blue}read.cross()}} used in 
the R package \texttt{qtl} \citep{rqtl}, and accepts the basic \texttt{R/qtl} CSV format. 
For more information about this function type \texttt{{\color{blue}?read.population}}. 
In this function the data file can be specified as an argument.\\

\texttt{> obesity.cross <- read.population("Obesity.Cross.csv")}\\

Alternatively, the filename can be left blank to choose a file through the file system. The following
code will bring up a window for choosing the desired file:\\

\texttt{> obesity.cross <- read.population()}\\

If the phenotypes are not specified in \texttt{\color{blue}{read.population()}} the default 
behavior is to read in all phenotypes. The functions \texttt{{\color{blue}delete.pheno()}} 
and \texttt{{\color{blue}select.pheno()}} can be used after reading in the data to narrow 
down the phenotypes for analysis. For \texttt{cape} analysis there must be 
between two and 12 phenotypes. At least two phenotypes must be used because the
fundamental concept of \texttt{cape} uses partial pleiotropy to refine models of epistasis. 
However, too many phenotypes may reduce power to detect interactions.\\ 

Upon being read in, the data are formatted and stored in an object. This object is referred
to as \texttt{data.obj} in the argument lists of most functions in \texttt{cape}. The main
functions in this analysis return \texttt{data.obj} with any results appended to it. The 
structure of this object can be viewed at any time by using the core R function
\texttt{{\color{blue}str()}}. For example, the data that we have just loaded is 
named \texttt{obesity.cross}, and its structure can be seen here:

<<>>=
str(obesity.cross)
@

Individual elements of this list can be accessed by name by using the following syntax:\\

\verb@data.obj$name@\\

When the cross is first loaded into \texttt{cape} it contains the following elements:

\begin{itemize}
\item \textbf{\$pheno} - The phenotype matrix in which individuals are in 
rows and phenotype values are in columns.
\item \textbf{\$geno} - The genotype matrix in which individuals are in rows
and genetic loci are in columns. Each cell of the genotype matrix contains the 
probability that individual $i$ has the reference allele at that position. By default, the 
reference allele is the allele that is alphabetically first if the original data file used letter 
designations for genotypes, or ``0" if the original data file used numeric genotype 
designations. In the example backcross, NON/Lt is the reference allele and the heterozygote
state is coded as the perturbation. Thus \texttt{cape} will model the effects of NZO/HILt variants.
\item \textbf{\$chromosome} - The chromosome on which each genetic locus
is situated
\item \textbf{\$marker.names} - The alphanumeric names of markers.
\item \textbf{\$marker.location} - The position of each locus on each chromosome
\end{itemize}

Most functions performed on \texttt{data.obj} will add elements to this list containing
the results of the analyses.

\section*{Manipulating the Data Object}
\addcontentsline{toc}{section}{Manipulating the Data Object}

There are a number of functions available for simple manipulation of the data object. For
example, after the data have been read in, phenotypes may be removed from the phenotype
matrix using \texttt{{\color{blue}delete.pheno()}}. To select specific phenotypes, use the
function \texttt{{\color{blue}select.pheno()}}.\\

<<>>=
obesity.cross <- select.pheno(obesity.cross, 
phenotypes = c("body_weight", "glucose", "insulin", "mom"))
@

In \texttt{cape} covariates are imported as phenotypes and must be reassigned as 
genetic markers. To code a variable, such as sex, experimental treatment, or another 
environmental variable, as a covariate, use the function \texttt{{\color{blue}create.covar()}} 
to transfer a variable from the phenotype matrix to the genotype matrix. Here we convert 
the factor maternal environment (``mom'') to a covariate.\\

<<>>=
obesity.cross <- create.covar(obesity.cross, "mom")
@

The data object can also be subset by chromosome (\texttt{{\color{blue}select.by.chr()}}) 
or by individual (\texttt{{\color{blue}select.by.ind()}}) before the analysis if desired. 
Individuals can be selected based either on phenotype or genotype value. For example,
to use only individuals with a plasma insulin level less than than 25 ng/mL:\\

\verb@> obesity.cross <- select.by.ind(obesity.cross, "pheno", "insulin < 25")@\\

Finally, it is also possible to remove all individuals from the cross with missing phenotype
data (\texttt{{\color{blue}remove.ind.with.missing.pheno()}}). 


\section*{Examining the Data}
\addcontentsline{toc}{section}{Looking at the Data}

Before proceeding with an analysis it is recommended that the data be examined by 
eye. The R package \texttt{qtl} has sophisticated plotting tools for examining genetic cross data,
which we do not try to duplicate here. It is straightforward enough, however, to examine 
phenotype distributions for normality, batch effects and other quality control issues using 
standard R functions. Using the function \texttt{{\color{blue}hist()}}, we can look as the 
distribution of each phenotype:


<<fig=TRUE,echo=FALSE,height=3>>=
layout(matrix(c(1:3), nrow = 1))
hist.data <- apply(matrix(c(1:dim(obesity.cross$pheno)[2]), nrow = 1), 2, function(x) 
hist(obesity.cross$pheno[,x], main = colnames(obesity.cross$pheno)[x], 
xlab = colnames(obesity.cross$pheno)[x]))
@

While body weight looks relatively normally distributed, glucose and insulin have
obviously non-normal distributions. 

Examination of the Q-Q plots of pairs of phenotypes, using the core R function 
\texttt{{\color{blue}qqplot()}}, can reveal phenotyping errors and other pathologies. Here we 
see a threshold and ceiling effect in the relationship between the distributions of insulin and 
the other two phenotypes.


<<fig=TRUE,echo=FALSE,height=3>>=
layout(matrix(c(1:3), nrow = 1))
i <- 1
while(i < dim(obesity.cross$pheno)[2]){
	j <- i + 1
	while(j <= dim(obesity.cross$pheno)[2]){
		qqplot(obesity.cross$pheno[,i], obesity.cross$pheno[,j], 
		xlab = colnames(obesity.cross$pheno)[i], ylab = 
		colnames(obesity.cross$pheno)[j], cex.lab = 1.5)
		j <- j + 1
	}
	i = i + 1
}
@

In general we recommend mean centering and normalizing all phenotypes before
proceeding with the analysis. Phenotype normalization can be achieved through log 
transformation, quantile normalization, or another method before the analysis. The 
function \texttt{{\color{blue}norm.pheno()}} uses quantile normalization to fit the 
phenotypes to a normal distribution. Briefly, this process sorts the values of the 
phenotype and replaces each with a corresponding value drawn from a normal
distribution with the same standard deviation and mean as the original distribution.
Mean centering subtracts the mean phenotype value from each phenotype value
yielding a distribution centered around 0.

<<>>=
obesity.cross <- norm.pheno(obesity.cross, mean.center = TRUE)
@

Plotting the histograms of the normalized data confirms the normalization. Only insulin still has
a ceiling effect, which cannot be removed by normalization because rank cannot be determined
among equal values.

<<fig=TRUE,echo=FALSE,height=3>>=
layout(matrix(c(1:3), nrow = 1))
hist.data <- apply(matrix(c(1:dim(obesity.cross$pheno)[2]), nrow = 1), 2, function(x) 
hist(obesity.cross$pheno[,x], main = colnames(obesity.cross$pheno)[x], 
xlab = colnames(obesity.cross$pheno)[x]))
@

We can also see that the Q-Q plots from before show the example phenotypes have been 
converted to the same distribution. The ceiling effect is still visible in the insulin measurement,
but this cannot be removed through normalization, and overall, the distribution is more similar
to those of the other phenotypes than before normalization. Knowing that this ceiling effect is 
present will be important in interpreting the results of the analysis. 


<<fig=TRUE,echo=FALSE,width=7,height=3>>=
layout(matrix(c(1:3), nrow = 1))
i <- 1
while(i < dim(obesity.cross$pheno)[2]){
	j <- i + 1
	while(j <= dim(obesity.cross$pheno)[2]){
		qqplot(obesity.cross$pheno[,i], obesity.cross$pheno[,j], 
		xlab = colnames(obesity.cross$pheno)[i], ylab = 
		colnames(obesity.cross$pheno)[j], cex.lab = 1.5)
		j <- j + 1
	}
	i = i + 1
}
@

At this point, if we decide to exclude insulin from the analysis because its distribution cannot
be normalized, we can simply remove it from the data object using \texttt{{\color{blue}delete.pheno()}}.\\

\verb@ > obesity.cross <- delete.pheno(obesity.cross, phenotypes = "insulin")@\\

\hspace{30 mm}

For the example here, however, we will include it in the analysis.


\section*{A Note on Phenotype Selection}
Phenotype selection is an important component of the \texttt{cape} analysis and should be given
considerable thought. This method relies on the selection of two or more phenotypes that have 
common genetic factors but are not identical across all individuals. Such phenotypes may describe
multiple aspects of a single complex trait, such as obesity or diabetes, and may encompass a 
combination of molecular phenotypes, such as plasma glucose levels, and phenotypes, such as
body weight, that are measured at the organismal level. The central assumption of this method is
that different genetic interactions found for a single gene pair in the context of different phenotypes 
represent multiple manifestations of a single underlying gene network. By measuring the 
interactions between genetic variants in different contexts we can gain a clearer picture of the 
network underlying statistical epistasis \citep{Carter:2012fd}.\\

The phenotypes in the \cite{Reifsnyder:2000is} data set are ideal for the \texttt{cape} analysis. 
They measure different aspects of the diabetes and obesity, two complex traits that are known to 
be related biologically and are highly correlated. The phenotypes themselves range from molecular 
phenotypes to organismal phenotypes. By examining the correlations between phenotypes, 
we can see that the phenotypes measured in this experiment are correlated, but not identical 
across all individuals. Ideally, phenotypes used in cape should have a Pearson correlation 
coefficient r between 0.4 and 0.8.

<<fig=TRUE,echo=FALSE,height=3>>=
layout(matrix(c(1:3), nrow = 1))
i <- 1
while(i < dim(obesity.cross$pheno)[2]){
	j <- i + 1
	while(j <= dim(obesity.cross$pheno)[2]){
		plot(obesity.cross$pheno[,i], obesity.cross$pheno[,j], 
		xlab = colnames(obesity.cross$pheno)[i], ylab = 
		colnames(obesity.cross$pheno)[j], main = paste("r =", 
		round(cor(obesity.cross$pheno[,i], obesity.cross$pheno[,j], 
		use = "complete.obs"), 2)), cex.lab = 1.2)
		j <- j + 1
	}
	i = i + 1
}
@

Note that all mice in the example backcross are male. For multisex populations, sex
is usually a covariate and correlation should be assessed for each sex separately.

\section*{Decomposing the Phenotypes}
\addcontentsline{toc}{section}{Decomposing the Phenotypes}

Although \texttt{cape} can find genetic variants associated with raw phenotypes, the analysis
was designed to work on composite traits called ``eigentraits." Eigentraits are calculated by 
factoring the matrix of phenotypes by singular value decomposition (SVD):\\

\begin{equation*}
Y = U \cdot V \cdot W^{T}
\end{equation*}

Where $Y$ is a matrix containing one column for each mean-centered, normalized phenotype
and one row for each individual. If $Y$ contains more individuals than phenotypes, the $U$ matrix
has the same dimensions as Y with each column containing one eigentrait. $V$ contains the
singular values, and $W^{T}$ contains the right singular vectors in rows. See \cite{Carter:2012fd}
for more details.\\

The SVD de-correlates the phenotypes concentrating phenotypic features into individual 
eigentraits. One benefit of this process is that variants that are weakly correlated to
several phenotypes due to common underlying processes may be strongly correlated 
to one of the eigentraits. This eigentrait captures the information of the underlying process, making 
strong main effects distributed between phenotypes easier to detect and identify as potential 
interaction loci and/or covariates. Thus, analysis of eigentraits is recommended over the 
analysis of raw traits.\\

To decompose the phenotypes to eigentraits use the function \texttt{{\color{blue}get.eigentraits()}}
If the phenotypes have not already been mean centered and normalized, this step can be performed
here. In this example, we do not carry out these steps, as they have already been performed. 

<<>>=
obesity.cross <- get.eigentraits(obesity.cross, scale.pheno = FALSE,
normalize.pheno = FALSE)
@

This function performs the SVD and returns the matrix of eigentraits in a new element of the list
called \texttt{ET}. The result of the decomposition can be viewed with the function 
\texttt{{\color{blue}plotSVD()}}.


<<fig=FALSE>>=
plotSVD(obesity.cross, orientation = "vertical")
@

%<<echo=FALSE, results=hide, print=FALSE>>==
%dev.copy2pdf(file = "svd.pdf")
%dev.off()
%@

\begin{center}
\includegraphics[width=3.5in]{svd.pdf}
\end{center}

In the example illustrated here, the first eigentrait captures more than 70\% of the variance
in the three phenotypes. This eigentrait describes the processes by which body weight, 
glucose levels, and insulin levels all vary together. The correlations between obesity and risk
factors for obesity, such as elevated insulin and fasting glucose levels are well known 
\citep{Permutt:2005ir,Das:2006bo,Haffner:2003kb}. The second eigentrait captures nearly 20\%
of the variance in the phenotypes. It captures the processes through which glucose and body
weight vary in opposite directions. This eigentrait may be important in distinguishing the 
genetic discordance between obesity and diabetes. While obesity is a strong risk factor for diabetes, 
not all those who are obese have diabetes, and not all those with diabetes are obese 
\citep{Permutt:2005ir,Burcelin:2002co}.\\

The third eigentrait is less interpretable biologically, as it describes the divergence of blood 
glucose and insulin levels. It may represent a genetic link between glucose and body weight 
that is non-insulin dependent. Because we are primarily interested 
in the connection between diabetes and insulin, we will use only the first two eigentraits for 
the analysis. In many cases in which more than two phenotypes are being analyzed, 
the first two or three eigentraits will capture the majority of the variance in the 
data and capture obvious features. Other eigentraits may capture noise or systematic 
bias in the data. Often the amount of total variance captured by such eigentraits is small, 
and they can be removed from the analysis.\\

Ultimately, there is no universal recipe for selecting which eigentraits should be included 
in the analysis, and the decision will be based on how the eigentraits contribute to the original 
phenotypes and how much variance in the data they capture.\\

To select eigentraits for the analysis, use the function \texttt{{\color{blue}select.eigentraits()}}.
Here we select the first two eigentraits.

<<>>=
obesity.cross <- select.eigentraits(obesity.cross, traits.which = c(1,2))
@


\section*{Single-Variant Scan}
\addcontentsline{toc}{section}{Single-Variant Scan}

Once the eigentraits for the analysis have been selected, the single-locus scan is run to
investigate how individual markers are associated with each eigentrait. Note that this
scan performs a linear regression at each marker. A more sophisticated single-locus analysis 
can be performed by \texttt{R/qtl} \citep{rqtl}. This single-marker scan in \texttt{cape} 
performs the following regression for each locus on each eigentrait: 

\begin{equation*}
U_{i}^{j} = \beta_{0}^{j} + x_{i}\beta^{j} + \epsilon_{i}^{j}
\end{equation*}

The index $i$ runs from 1 to the number of individuals, and $j$ runs from 1 to the number
of eigentraits or phenotypes. $x_{i}$ is the probability of the presence of the reference allele for 
individual $i$ at locus $j$. The purpose of this scan is two-fold:

\begin{itemize}
\item In large data sets the number of possible variant pairs may be too large to test exhaustively.
The single-variant scan can be used as a filtering step to choose variants that will be included in
the pair scan.
\item Large main effects can obscure interactions. In sufficiently powered studies, conditioning
 on the large QTL can aid in the 
discovery of interactions, variants with large main-effects can be used as covariates in the pair scan.
\end{itemize}

The function \texttt{{\color{blue}singlescan()}} runs the single-marker scan and has several 
arguments associated with selecting markers to include in the pairwise scan as well as for 
selecting markers to use as covariates in the pairwise scan.

<<>>=
obesity.cross <- singlescan(obesity.cross, n.perm = 100, covar = "mom", 
scan.what = "eigentraits", alpha = c(0.01, 0.05), verbose = FALSE)
@

The argument \texttt{alpha} takes a vector of alpha values for which significance thresholds
will be calculated. These values are used primarily for plotting purposes. 

%Markers can be selected for inclusion in the pairwise scan based on their having significant
%main effects as determined through permutation testing. The significance threshold is set with 
%the argument \texttt{alpha.for.pairs}. After the single-locus scan, if filtering markers by 
%significance is desired, the function \texttt{{\color{blue}select.markers.for.pairscan}} selects
%only those markers exceeding this significance cutoff for inclusion in the pairscan. If no filtering
%is desired, this threshold is ignored.\\

%The function \texttt{{\color{blue}singlescan()}} can also automatically select markers to
%include as covariates in the pairwise scan. Like the marker inclusion selection, marker 
%covariate selection is also based on marker effect size, and the alpha for this threshold 
%is set by the argument \texttt{alpha.for.covar}. Automatic covariate selection is recommended 
%in populations in which genetic markers are independent of one another, but not in populations 
%in which markers are linked (see \underline{\nameref{sec:Covariate Selection}}). To use 
%automatic covariate selection, set the argument \texttt{auto.covar.selection} to 
%\texttt{TRUE}.\\

Covariates to be used in the single-marker scan should be specified through the argument
\texttt{covar}. These covariates should have already been moved to the marker matrix
by \texttt{{\color{blue}create.covar()}}. An additional argument \texttt{scan.what} determines 
whether ``eigentraits" or ``raw.traits" will be used as the dependent variable in the regression.\\

The single-marker scan currently does not support markers on sex chromosomes. Because the X
chromosome is hemizygous in males, sex differences in phenotype can lead to false associations,
and markers on this chromosome require special consideration \citep{Broman:2006fb}. Before the
single-marker scan is performed, any markers on the X and Y chromosomes are removed from the 
data object. The results of the single variant scan can be visualized with the function
\texttt{{\color{blue}plotSinglescan()}}.

<<fig=TRUE, width=6>>=
plotSinglescan(obesity.cross, mark.chr = TRUE, mark.covar = FALSE)
@

In this figure the t-statistic ($\beta/\sigma$) of each marker is plotted as a vertical line.
Results for both eigentraits are shown here as ET1 and ET2, and chromosome numbers are
written along the x axis. If filtering of markers for the pair-wise scan is desired, an effect size
cutoff can be specified using either \texttt{{\color{blue}set.pairscan.thresh()}} or in 
\texttt{{\color{blue}select.markers.for.pairscan()}}. If a specific number of markers is desired, 
\texttt{{\color{blue}select.markers.for.pairscan()}} can determine the threshold at which this number
of markers is included in the pair scan. Setting a threshold for inclusion is useful for large data 
sets for which it is impractical to test all marker pairs. Here we use all markers in the pairwise scan. 
In this example we will not filter the markers.\\

%The higher dashed line represents the alpha value defined by the argument \texttt{alpha.for.covar},
%which here was set to the default value of 0.01. This threshold can be used for automatic covariate
%selection. If the argument \texttt{automatic.covar.selection} is set to \texttt{TRUE} all markers 
%exceeding the \texttt{alpha.for.covar} threshold will automatically be designated as covariates. This 
%automated selection is generally only recommended for populations in which variants 
%segregate independently. For more details on covariate selection, see 
%\underline{\nameref{sec:Covariate Selection}}.


\section*{Covariate Selection}
\addcontentsline{toc}{section}{Covariate Selection}
\label{sec:Covariate Selection}
As mentioned previously, conditioning on large main effects may aid in the discovery of
interactions. Therefore \texttt{cape} allows specification of genetic markers as covariates.
%which can be done either manually or automatically. Automatic covariate selection is 
%particularly useful in populations in which variants segregate independently. For example, 
%\emph{Drosophila} have virtually no linkage between markers \citep{Mackay:2012fd}; 
%therefore most markers can be assumed to assort independently. Similarly, mutations 
%to genes on different chromosomes in a yeast population segregate independently. 
%In these and other similar cases, automatic covariate selection is appropriate because 
%markers achieving the significance threshold can be assumed to be independent of each
%other and will not reduce the significance of other markers due to linkage if used as covariates.\\

This specification is recommended if there are markers with particularly large effect sizes
that may obfuscate weak interactions. In an organism in which variants segregate relatively
independently, such as in \emph{Drosophila} \citep{Mackay:2012fd}, it may be useful to
select individual markers with very high effect sizes as covariates. In a population with more
correlation between adjacent markers, more care should be given to this process. In a mouse 
F2 intercross or backcross, on the other hand, markers on a single chromosome tend to be 
linked, and groups of markers may be significantly associated with a phenotype due to linkage 
with a causative locus. In this case, markers from a single linkage block carry partially redundant 
information and, if treated as independent covariates, will obfuscate interactions between makers 
in the linkage block and other truly independent markers. Thus, if markers in the study population 
are linked, it is recommended that any marker covariates be set manually and that only one 
covariate should be selected from each linkage block.\\

Setting covariates, as well as the filtering of markers for the pair-wise scan can be done with 
the function \texttt{{\color{blue}select.markers.for.pairscan()}}. When calling this function, if
thresholding the markers based on effect size, set \texttt{use.pairs.threshold} to TRUE, and
\texttt{pairscan.thresh} to the desired standardized effect size threshold. Alternatively, this
function can be used to select a threshold based on a desired number of markers to test
in the pair-wise scan. The target number of markers desired is set in \texttt{num.markers}.
The function then searches through multiple thresholds to find a linearly independent 
matrix holding the number of desired markers. It starts the search at the threshold set in
\texttt{start.thresh} and stops when it finds a threshold resulting in \texttt{num.markers} 
plus or minus the \texttt{tolerance}. This process is slow, and it is recommended that the 
starting threshold be set as close to the estimated correct threshold as possible and run
only once to find the desired threshold. After the first run, the threshold can be set manually.

In addition to selecting markers based on an effect size threshold, 
\texttt{{\color{blue}select.markers.for.pairscan()}} can also set specific markers as covariates
through the argument \texttt{specific.markers}. This argument takes a vector of strings and
sets each one as a covariate. 

<<>>=
obesity.cross <- select.markers.for.pairscan(obesity.cross, 
use.pairs.threshold = FALSE)
@


%To use automatic covariate selection, set \texttt{auto.covar.selection} to \texttt{TRUE} when
%running \texttt{{\color{blue}}singlescan}. All markers above the \texttt{alpha.for.covar} threshold 
%will be automatically included as covariates in the pair scan. If this option is chosen, is is
%still possible to adjust the covariate threshold after the single-marker scan to include
%more or fewer covariates as desired. In general we recommend that only the 
%markers with the very strongest effects are treated as covariates to avoid overfitting the 
%model. The function \texttt{{\color{blue}get.covar()}} sets the covariate threshold in the 
%data object and recalculates which markers are to be used as covariates in the pair scan.\\
%
%
%\verb@ > obesity.cross <- get.covar(obesity.cross, covar.thresh = 4.5)@\\
%
%If selecting covariates manually, set \texttt{auto.covar.selection} to \texttt{FALSE} when
%running \texttt{{\color{blue}}singlescan}. The plots from the single-marker scan can
%then be examined to determine which markers should be used as covariates. The 
%function \texttt{{\color{blue}set.covar()}} can then be used to set covariates individually.\\
%
%Here, because of the large effects of maternal environment on ET1, we make this marker 
%a covariate. The function \texttt{{\color{blue}set.covar()}} optionally includes a plot to verify 
%the covariate selection. Covariates are marked in red.
%\\
%
%<<echo=TRUE,fig=TRUE>>=
%obesity.cross <- set.covar(obesity.cross, pheno = "ET1",  
%markers = c("mom"), is.covar = TRUE, plot.covar = TRUE)
%@



\section*{Pairwise Scan}
\addcontentsline{toc}{section}{Pairwise Scan}

The purpose of the pairwise scan is to find interactions, or epistasis, between variants. 
The epistatic models are then combined across phenotypes or eigentraits 
to infer a parsimonious network that takes data from all eigentraits into account.\\

To find epistatic interactions \texttt{{\color{blue}pairscan()}} tests the following model for each 
variant 1 and 2:

\begin{equation*}
U_{i}^{j} = \beta_{0}^{j} + \underbrace{\sum_{c=1}^{2}x_{c,i}\beta_{c}^{j}}_{\mathrm{covariates}} + 
\underbrace{x_{1,i}\beta_{1}^{j} + x_{2,i}\beta_{2}^{j}}_{\mathrm{main\;effects}} + \underbrace{x_{1,i}x_{2,i} 
\beta_{12}^{j}}_{\mathrm{interaction}} +  \epsilon_{i}^{j}
\end{equation*}

The terms in this equation are the same as those in the equation for the single-variant scan
except for the addition of the term for the interaction between the two variants being tested.
This additional term brings further complications to the model, and restricts which 
markers can be tested. Because many markers, including covariates, may be included in 
the model, we need to be careful about including only markers that are linearly independent 
of each other. Linear dependence between markers occurs when two markers are fully 
linked and therefore perfectly correlated. The two markers provide the same genetic information, 
and one can be discarded without loss of information. Before running the pair scan, it is important 
that we reduce the genetic matrix to only markers that are linearly independent of one another. 
This step is performed by the function \texttt{{\color{blue}select.markers.for.pairscan()}}. It
calculates the correlation between all pairs of markers. If any are found to be perfectly correlated, 
the first marker is removed.\\

As mentioned above, \texttt{{\color{blue}select.markers.for.pairscan()}} also optionally filters 
markers for inclusion in the pair scan by standardized effect size. This is useful in large 
crosses in which it may be impossible to test all possible pairs of markers in a reasonable 
amount of time.

\begin{verbatim}
obesity.cross <- select.markers.for.pairscan(obesity.cross, 
use.pairs.threshold = FALSE)
\end{verbatim}

The number of markers removed is printed to the screen, and the names of these markers
are written to the file \texttt{markers.removed.txt} in the current working directory. This filtering
step adds two elements to the data object. One is a filtered genotype matrix called 
\texttt{geno.for.pairscan}. The other, \texttt{covar.for.pairscan}, indicates which markers are 
to be used as covariates in the pair scan.\\

The results of the filtering can be visualized with
\texttt{{\color{blue}ploSinglscan()}}. This visualization can be helpful in seeing
the locations of discarded markers, especially if many have been removed. To
indicate which markers have been selected for the pair scan, set 
\texttt{show.selected.markers} to \texttt{TRUE}. To indicate which markers have 
been rejected from the pair scan, set \texttt{show.rejected.markers} to \texttt{TRUE}.

<<fig=TRUE>>=
plotSinglescan(obesity.cross, mark.chr = TRUE, show.rejected.markers = TRUE, 
standardized = TRUE)
@


After ensuring that all markers are linearly independent and thresholded satisfactorily, the
pair scan can be run using \texttt{{\color{blue}pairscan()}}.


\begin{verbatim}
obesity.cross <- pairscan(obesity.cross, scan.what = "eigentraits", 
n.perm = 100, min.per.genotype = 6, verbose = FALSE)
\end{verbatim}

%<<echo=FALSE, results=hide, print=FALSE>>=
%obesity.cross <- pairscan(obesity.cross, scan.what = "eigentraits", 
%n.perm = 10, min.per.genotype = 6, verbose = FALSE)
%@

The arguments here are familiar from \texttt{{\color{blue}singlescan()}}, with the exception
of two additional arguments. \texttt{min.per.genotype} and \texttt{max.pair.cor} are thresholds 
that prevent highly correlated markers from being tested in pairs. Only one of these arguments 
can be set depending on whether the genotypes in the data set are continuous or discrete.
If the genotypes are discrete, \texttt{min.per.genotype} can be set to prevent empty cells in 
genotype matrices. In an intercross, there are three possible genotypes at each marker, 
giving a total of nine possible genotypes for the pair. For example, for two markers each 
with the genotypes AA, AB, and BB, the pairwise genotypes are the following:\\

\begin{center}
\begin{tabular}{ r|c|c|c| }
\multicolumn{1}{r}{}
 &  \multicolumn{1}{c}{AA}
 & \multicolumn{1}{c}{AB}
 & \multicolumn{1}{c}{BB}\\
\cline{2-4}
AA & & & \\
\cline{2-4}
AB & & & \\
\cline{2-4}
BB & & & \\
\cline{2-4}
\end{tabular}
\end{center}
\hspace{30 mm}

In a backcross each marker only has two possible genotypes, AA and AB, yielding four 
possible pairwise genotypes. In both cases, insufficient representative individuals of single 
genotypes indicates that the two markers in question are linked. This linkage leads to false 
associations in both the effects and permutations. False associations in permutations results 
in a heavy-tailed null distribution, which artificially inflates the threshold of significance and 
reduces power to find true interactions. Because of this effect, linked marker pairs are not tested. 
The threshold used to determine insufficient recombination between markers is given by 
\texttt{min.per.genotype}. The default behavior is to reject any marker pair for which there are
fewer than six individuals in any of the genotype cells. This value can be adjusted, but caution 
should be used in interpreting the results if \texttt{min.per.genotype} is very low or 0.\\

Alternatively, if genotypes are coded continuously, the argument \texttt{max.pair.cor} can be 
set to establish an upper threshold on the Pearson correlation between marker pairs. As in 
the case with discrete genotypes, testing highly correlated markers against each other can
lead to false associations. Either \texttt{max.pair.cor} or \texttt{min.per.genotype}, but not both,
must be set.\\

The \texttt{n.perm} argument in \texttt{{\color{blue}pairscan()}} sets the number of 
permutations performed per marker pair. To generate a null distribution for the 
pair-wise tests, CAPE performs permutation tests in the following way: First a singlescan
is run on the permuted phenotypes. Then the N markers with the largest effect sizes 
are selected, where N is the same number of markers being tested in the pairwise scan.
This ensures that if there are very few markers being tested in the pairwise scan, we do not
generate a biased null distribution by only performing permutation tests repeatedly on 
these markers. A pairwise scan is then performed on these top N markers using the permuted
phenotype. To calculate empirical p values, the permutations are combined across all markers 
to result in one large null distribution. Thus, large null distributions can be achieved with 
relatively few permutations per pair. This is useful, as even with relatively sparse genotyping 
the number of tests performed can be large and take many hours to perform.\\

The results of the pair scan can be plotted with \texttt{{\color{blue}plotPairscan()}}.
This plot shows the resulting $\beta/\sigma$ for each pair of markers. Gray and 
white bars show the boundaries of the chromosomes.

%<<fig=FALSE, results=hide, print=FALSE>>=
%plotPairscan(obesity.cross, phenotype = c("ET1", "ET2"), standardized = FALSE, 
%pdf.label = "Pairscan.Regression.pdf")
%@

\begin{verbatim}
plotPairscan(obesity.cross, phenotype = c("ET1", "ET2"))
\end{verbatim}
\includegraphics{pair_regression.pdf}


\section*{Combined Analysis for Detection of Interactions}
\addcontentsline{toc}{section}{Combined Analysis for Detection of Interactions}

From the pair scan, each pair of markers 1 and 2 receives a set of $\beta$ coefficients
describing the main effect of each marker on each eigentrait $j$ ($\beta^j_1$ and 
$\beta^j_2$) as well as the interaction effect of both markers on each eigentrait 
($\beta^j_{12}$) (See figure below). The central idea of \texttt{cape} is that these 
coefficients can be combined across eigentraits and reparameterized to calculate 
how each pair of markers influences each other directly and independently of eigentrait.

\begin{center}
\label{fig:reparam}
\includegraphics[width=4in]{reparam.pdf}
\end{center}

The first step in this reparameterization is to define two new parameters ($\delta_1$ 
and $\delta_2$) in terms of the interaction coefficients. $\delta_1$ can be thought of 
as the additional genetic activity of marker 1 when marker 2 is present. Together the
$\delta$ terms capture the interaction term, and are interpreted as the extent to which 
each marker influences the effect of the other on downstream phenotypes. For example, 
a negative $\delta_2$ indicates that the presence of marker 2 represses the effect of 
marker 1 on the phenotypes or eigentraits. The $\delta$ terms are related to the main 
effects and interaction effects as follows:

\begin{equation}
\label{eqn:delta_mat_def}
\begin{bmatrix}
\beta^1_1 & \beta^1_2\\
\beta^2_1 & \beta^2_2\\
\vdots & \vdots
\end{bmatrix}
\cdot
\begin{bmatrix}
\delta_1\\
\delta_2\\
\end{bmatrix}
=
\begin{bmatrix}
\beta^1_{12}\\
\beta^2_{12}\\
\vdots
\end{bmatrix}
\end{equation}

In multiplying out this equation, it can be seen how the $\delta$ terms influence
each main effect term to give rise to the interaction terms independent of phenotype. 

\begin{equation}
\label{eqn:delta_def}
\beta^j_1\delta_1 + \beta^j_2\delta_2 = \beta^j_{12}
\end{equation}

If $\delta_1 = 0$ and $\delta_2 = 0$, there are no addition effects exerted by the 
markers when both are present. Substitution into the equations above shows that 
the interaction terms $\beta^j_{12}$ are $0$ and thus the interaction terms have 
no effect on the phenotype.\\
 
Alternatively, consider the situation when $\delta_1 = 1$ and $\delta_2 = 0$. The positive
$\delta_1$ indicates that marker 1 should exert an additional effect when marker 2 is 
present. This can be seen again through substitution into equation \ref{eqn:delta_def}:

\begin{equation*}
\beta_j^1 = \beta_{12}^j
\end{equation*}

These non-zero terms show that there is an interaction effect between marker 1 and 
marker 2. The positive $\delta_1$ indicates that this interaction is driven through
an enhanced effect of marker 1 in the presence of marker 2.\\

The $\delta$s are calculated by solving for equation \ref{eqn:delta_mat_def} using 
matrix inversion:

\begin{equation*}
\begin{bmatrix}
\delta_1\\
\delta_2\\
\end{bmatrix}
=
\begin{bmatrix}
\beta^1_1 & \beta^1_2\\
\beta^2_1 & \beta^2_2\\
\vdots & \vdots
\end{bmatrix}^{-1}
\cdot
\begin{bmatrix}
\beta^1_{12}\\
\beta^2_{12}\\
\vdots
\end{bmatrix}
\end{equation*}

This inversion is exact for two eigentraits, and \texttt{cape} implements pseudo-inversion 
for up to 12 eigentraits.\\

The $\delta$ terms are then translated into directed variables defining the 
marker-to-marker influences $m_{12}$ and $m_{21}$. Whereas $\delta_2$ described the
change in activity of marker 2 in the presence of marker 1, $m_{12}$ can be thought 
of as the direct influence of marker 2 on marker 1, with negative values indicating repression
and positive values indicating enhancement. The terms $m_{12}$ and $m_{21}$ are
self-consistent and defined in terms of $\delta_1$ and $\delta_2$:

\begin{equation*}
\delta_1 = m_{12}(1 + \delta_2),\;\delta_2 = m_{21}(1 + \delta_1)
\end{equation*}
 
Rearranging these equations yields the solutions:
 
\begin{equation*}
m_{12} = \frac{\delta_1}{1 + \delta_2},\;m_{21} = \frac{\delta_2}{1 + \delta_1}.
\end{equation*}
 
These directed influence variables provide a map of how each marker influences each
other marker independent of phenotype. The significance of these influences is determined
through standard error analysis on the regression parameters \citep{bevington1994data,Carter:2012fd}.
This step is particularly important as matrix inversion can lead to large values but larger
standard errors, yielding insignificant results. As an example, the variance of $m_{12}$
is calculated by differentiating with respect to all model parameters:

\begin{equation*}
\sigma^{2}_{m_{12}} \cong \sum_{ij} \sigma^2_{\beta_{i}^{j}} \Bigg(\frac{\partial m_{12}}
{\partial \beta_{i}^{j}}\Bigg)^2 + 2 \sum_{i<k, j < l} \sigma^2_{\beta^{j}_{i}\beta^{l}_{k}} \Bigg(\frac{\partial m_{12}}
{\partial \beta^{j}_{i}} \Bigg) \Bigg(\frac{\partial m_{12}}{\partial \beta^{l}_{k}} \Bigg)
 \end{equation*}
 
In this equation, the indices $i$ and $k$ run over regression parameters and $j$ and $l$
run from 1 to the number of traits. The calculations of the $\delta$ and the $m$ terms, as well 
as the error propagation are performed by the function \texttt{{\color{blue}error.prop()}}. 
 
%<<>>=
%obesity.cross <- error.prop(obesity.cross, perm = FALSE, verbose = FALSE)
%@
 
\begin{verbatim}
obesity.cross <- error.prop(obesity.cross, perm = FALSE, verbose = FALSE)
\end{verbatim}

 This function is applied to both the results from the pairwise scan, as well as the permutations 
 of the pairwise scan for later calculation of p values. To apply the calculations to the pairwise
 scan permutations, set \texttt{perm = TRUE}.
 
%<<>>=
%obesity.cross <- error.prop(obesity.cross, perm = TRUE, verbose = FALSE)
%@

 \begin{verbatim}
obesity.cross <- error.prop(obesity.cross, perm = TRUE, verbose = FALSE)
\end{verbatim}
 
For large scans it might be desirable to observe the progress of the calculations. This can 
be done by setting \texttt{verbose = TRUE}. After these calculations have been performed
the results of the permutation testing can be used to calculate empirical p values for each
of the variant-to-variant effects.


%<<>>=
%obesity.cross <- calc.p(obesity.cross, pval.correction = "fdr")
%@

\begin{verbatim}
obesity.cross <- calc.p(obesity.cross, pval.correction = "fdr")
\end{verbatim}


This function also adjusts the empirical p values for multiple testing. The default correction 
is Holm's stepdown procedure \citep{holm1979simple}. Two other methods, false discovery rate 
(FDR) \citep{fdr} and local false discovery rate (lFDR) \citep{lfdr} are also available. The latter methods
of correction for multiple testing are less stringent than the Holm's step-down procedure. The function 
\texttt{{\color{blue}calc.p()}} adds two elements to the data object. The purpose and returned results of 
each function is summarized below. For tables of all functions, see the 
\underline{\nameref{sec:Tables of Functions}} section at end of this document.\\

\begin{table}[h]\footnotesize
\caption*{\small{\textbf{Functions for Parsing Pair Scan Results}}}
\begin{tabular}{l P{3.5cm} P{3.5cm} P{2.5cm} P{3cm}}
\scriptsize{\textbf{Function}} & \scriptsize{\textbf{Behavior}} & \scriptsize{\textbf{Result Label}} & \scriptsize{\textbf{Structure}} & \scriptsize{\textbf{Columns}}\\
\hline \hline
\nohyphens
error.prop & propagates errors of coefficients calculated in the pair scan & var.to.var.influences or var.to.var.influences.perm & matrix & marker1, marker2, m12, m21, m12~$\sigma$, m21~$\sigma$\\
\hline
calc.p & calculates empirical and corrected p~values or fdr of interactions using permutations & empirical.p & list of 2 matrices (m12, m21)  & marker1, marker2, source, target, influence coefficient, influence $\sigma$, empirical p~value, adjusted p~value or fdr\\
\hline
\end{tabular}
\end{table}


\normalsize

The final step in calculating the network of directed influences is to translate the effect of each 
marker on the eigentraits to effects on the original phenotypes. The effects of variants on eigentraits
are defined in terms of the interaction coefficients $m_{12}$ and $m_{21}$ as well as the main effect 
of each variant $\beta_{i}^{j}$. To translate these effects to be in terms of the original phenotypes, 
the coefficient matrices are multiplied by the singular value matrices $V \cdot W^{T}$. With two 
phenotypes and two eigentraits this conversion results in no loss of information. It should also be
noted that the translation back to phenotype space does not affect the variant-to-variant influences.\\

The translation to phenotype space is performed by the function \texttt{{\color{blue}direct.influence()}}.

%<<>>=
%obesity.cross <- direct.influence(obesity.cross, pval.correction = "fdr")
%@

\begin{verbatim}
obesity.cross <- direct.influence(obesity.cross, pval.correction = "fdr")
\end{verbatim}


In these calculations each marker is assigned multiple direct influences on each phenotype
dependent on the marker it was paired with. To reduce the results to a single direct influence
of each marker on each phenotype, \texttt{{\color{blue}direct.influence()}} selects the maximum
influence observed for each marker. These effects are stored in \texttt{max.var.to.pheno.influence}.
All effects are accessible in the element \texttt{variant.to.phenotype.influences}. To save memory, 
data from the permutations are not saved by default; however, if \texttt{save.permutations} is set to
\texttt{TRUE} all data are saved in the object \texttt{permutation.data.RData}. For details on the
elements added to \texttt{obesity.cross} and those saved in \texttt{permutation.data.RData}, see 
the \underline{\nameref{sec:Tables of Functions}} section at end of this document.\\

The function \texttt{{\color{blue}direct.influence()}} also performs correction for multiple testing
for each of the influences. Holm step-down correction \citep{holm1979simple}, FDR \cite{fdr},
and local FDR \citep{lfdr} are all available as correction methods.\\


The marker-to-marker and marker-to-phenotype influences can be plotted with the function 
\texttt{{\color{blue}plotVariantInfluences()}}. This function plots the adjacency 
matrix of the final network. It shows all significant influences between variants and between
variants and phenotypes.


%<<>>=
%plotVariantInfluences(obesity.cross, p.or.q = 0.01, 
%all.markers = FALSE, standardize = FALSE, not.tested.col = "lightgray")
%@

%<<echo=FALSE, results=hide, print=FALSE>>==
%dev.copy2pdf(file = "variant_influences.pdf")
%dev.off()
%@

\begin{verbatim}
plotVariantInfluences(obesity.cross, p.or.q = 0.01, 
all.markers = FALSE, standardize = TRUE, not.tested.col = "lightgray", 
scale.effects = "log10",  pheno.width = 8)
\end{verbatim}


\begin{center}
\label{fig:varinf}
\includegraphics[width=6in]{variant_influences.pdf}
\end{center}
 
 
There are several arguments in this function to note. The argument \texttt{p.or.q} takes in
the p, q, or local fdr value at which interactions are considered significant. Only significant 
interactions are plotted. Non-significant interactions are colored white. It should be noted that 
some of the white boxes have gray dots in them. These dots indicate that these pairs were not 
tested for interactions because they were filtered out due to insufficient representation of individual 
phenotypes. This usually occurs because the markers are linked and there is low recombination 
between them. These pairs are marked with dots to distinguish them from pairs that were tested 
but did not have significant interactions between them. The color of the dots can be changed 
with the argument \texttt{not.tested.col}. The default is \texttt{lightgray}. This can be changed to 
\texttt{FALSE} or \texttt{NA} if no distinction between not-tested and not-significant is desired. 
Gray and white bars along the margins of the plot indicate the boundaries of the chromosomes.\\

It should also be noted that not all of the original markers from the single-locus scan are represented in
the final figure. This figure is showing only those markers that were included in the pair scan. The figure
can be changed to show all markers by setting \texttt{all.markers} to TRUE. However, with large numbers 
of markers, especially if not all markers are included in the pair scan, the figure can get quite crowded if all 
markers are included.\\

The argument \texttt{scale.effects} allows application of a function, either $log_{10}$ or square root
to the values in the adjacency matrix. These transformations can help increase the contrast between
markers. And finally, the argument \texttt{pheno.width} specifies how wide the main effects should
be relative to the interaction effects for better visualization.\\


A table of significant influences can be written to a file using the function\\
\texttt{{\color{blue}writeVariantInfluences()}}.

%<<>>=
%writeVariantInfluences(obesity.cross, p.or.q = 0.05, 
%filename =  "Significant.Influences.txt")
%@

\begin{verbatim}
writeVariantInfluences(obesity.cross, p.or.q = 0.05, 
filename =  "Significant.Influences.txt")
\end{verbatim}


\section*{Interpretation of Results}
\addcontentsline{toc}{section}{Interpretation of Results}

The table of variant-to-variant and variant-to-phenotype influences is the primary output
of \texttt{cape}. The function for plotting results, \texttt{{\color{blue}plotVariantInfluences()}},
plots an asymmetric adjacency matrix, which is useful for identifying patterns in
variant interactions.\\

Another useful visualization of the network can be plotted using \texttt{{\color{blue}plotNetwork()}}.
This function plots the chromosomes in a circle and shows interactions as arrows between
regions on chromosomes (see figures below). Main effects are plotted in the circles around the 
chromosomes. Before plotting the network in this format, the function \texttt{{\color{blue}get.network()}} 
must be be run to create a network object. The function \texttt{{\color{blue}get.network()}} network offers 
the option of condensing markers based on their linkage. The minimum Pearson $r$ 
threshold separating two markers is provided in the argument \texttt{r.thresh}.
The argument \texttt{p.or.q} determines the p, q, or local FDR value at which marker influences
are determined significant.


%<<>>==
%obesity.cross <- get.network(obesity.cross, p.or.q = 0.01, 
%collapse.linked.markers = TRUE, r2.thresh = 0.8, standardize = FALSE)
%plotNetwork(obesity.cross, collapsed.net = TRUE)
%@
%
%<<echo=FALSE, results=hide, print=FALSE>>==
%dev.copy2pdf("Network_Collapsed_r2.pdf")
%dev.off()
%@

\begin{verbatim}
obesity.cross <- get.network(obesity.cross, p.or.q = 0.01, 
collapse.linked.markers = TRUE, r.thresh = 0.5)
plotNetwork(obesity.cross, collapsed.net = TRUE)
\end{verbatim}


\begin{figure}[!ht]
\begin{center}
\includegraphics[width=5in]{Network_Collapsed.pdf}
\end{center}
\end{figure}

\pagebreak

If no collapse by linkage is desired, the argument \texttt{collapse.linked.markers} can be set
to \texttt{FALSE} and \texttt{r.thresh} will be ignored.


%<<>>==
%obesity.cross <- get.network(obesity.cross, p.or.q = 0.01, 
%collapse.linked.markers = FALSE, standardize = FALSE)
%plotNetwork(obesity.cross, collapsed.net = FALSE)
%@
%
%<<echo=FALSE, results=hide, print=FALSE>>==
%dev.copy2pdf("Network_Collapsed_full.pdf")
%dev.off()
%@

\begin{verbatim}
obesity.cross <- get.network(obesity.cross, p.or.q = 0.01, 
collapse.linked.markers = FALSE)
plotNetwork(obesity.cross, collapsed.net = FALSE)
\end{verbatim}


\begin{figure}[!ht]
\begin{center}
\includegraphics[width=5in]{Network_Full.pdf}
\end{center}
\end{figure}


Both the network figure and the adjacency matrix show direct influences of 
markers on the phenotypes as well as interactions between markers. As 
expected, NZO variants on multiple chromosomes show positive effects on 
plasma insulin and glucose levels as well as on body weight. Maternal 
obesity was found to suppress the effects of markers on Chromosome 15, 
and a marker on Chromosome 2 suppressed the effects of a marker on 
Chromosome 18. In the latter interaction, neither marker had an individual 
main effect, although the Chromosome 18 marker has a marginally significant 
effect (q = 0.018) on body weight that is completely suppressed by the 
Chromosome 1 marker.
%In R/cape markers are assigned significant main effects only
%if they alter phenotype independent of pair context. This reduces the incidence of false 
%positive main effects that depend on genetic background. While the Chromosome 2 
%and Chromosome 18 markers do not independently influence any of the phenotypes, 
%it can be seen from the interaction diagrams (Figure \ref{fig:interaction}) that the marker 
%at Chromosome 18 has a positive influence on body weight when the NON allele is 
%present on chromosome 2. However, the presence of the NZO allele on Chromosome 
%2 completely suppresses this positive effect.
NZO variants in a region on Chromosome 1 increased all three phenotypes individually 
 and were the hub of a small epistatic network suppressing the effects of markers on 
 Chromosomes 10 and 12.\\

The interaction between Chromosomes 1 and 12 is illustrative of the R/cape
 strategy. Although both loci have a significant effect on body weight, their joint effect is less
 than additive (See figure below). This finding suggests one locus might be 
 suppressing the effect of the other, but is ambiguous about the potential direction of 
 that suppression. In contrast, only the Chromosome 1 locus has an effect on 
 insulin and this effect is independent of the Chromosome 18 locus. This second 
 phenotype provides R/cape the information necessary to infer the directionality 
 of the interaction from Chromosome 1 to 12, since a reversed interaction would
 imply epistasis for insulin along with body weight, which was not observed.
 R/cape thus provides a more stringent hypothesis for gene candidates in the two
 loci through the constraint of directional genetic effects.\\

\begin{figure}[!ht]
\begin{center}
\includegraphics[width=6in]{interaction_fig.pdf}
\end{center}
%\caption{{\bf R/cape-derived interaction between Chromosome 1 and Chromosome 12.}
%(A) Effects of each genotype combination on each phenotype. The Chromosome 1 marker positively
%affects both insulin and weight, whereas the Chromosome 12 marker only has a significant effect
%on body weight. Note that plotted data is not conditioned on maternal obesity, which was a 
%covariate in the analysis.
%(B) Isolated interaction network for these loci, taken from Figure  \ref{fig:figresults}D, showing
%interaction between the loci (red edge between the circular nodes labeled by chromosome) 
%and their effects on the
%phenotypes (square nodes with labels "I" for insulin and "W" for weight).
%}
\label{fig:interaction}
\end{figure}


%The covariate maternal obesity (``mom'') is included in the network and was 
%found to suppress the effects of markers on chromosome 15. An interaction between a 
%chromosome 2 marker and a chromosome 18 marker highlight the ability of this method to 
%find indirect phenotypic modifiers. The chromosome 2 marker was found to suppress the effects 
%of a marker on chromosome 18 on body weight. However, the marker on chromosome 2 had no 
%significant main effect on body weight independent of the marker on chromosome 18. Its effects 
%on body weight are exerted through its interaction with the marker on chromosome 18. The model 
%is therefore that an NZO allele at the chromosome 2 locus indirectly affects body weight by 
%diminishing the weight-increasing effect of an NZO allele on chromosome 18.\\
%
%Loci on chromosome 1 suppress the effects of loci on chromosomes 3, 7, and 10 forming
%a small epistatic network with chromosome 1 at the center. NZO alleles on chromosomes 
%1, 3 and 10 were found individually to have positive effects on each of the phenotypes, 
%but \texttt{cape} determined that these effects were suppressed in animals with NZO alleles 
%at each pair of interacting locations. This pattern of interactions suggests redundancy 
%in the effects of multiple loci. While an NZO allele on chromosome 3 is hypothesized 
%to increase plasma glucose levels, it is hypothesized that an additional NZO allele on 
%chromosome 1 will not additively increase glucose levels further. By assigning directionality 
%and a quantitative magnitude to the genetic interaction, \texttt{cape} provides more specific 
%hypotheses of genetic interaction than standard epistasis analysis and uncovers indirect 
%genetic influences.\\


These findings illustrate how \texttt{cape} is designed to find interactions that 
simultaneously model all phenotypes under the assumption that interactions 
between variants across multiple contexts represent a single underlying 
interaction network. Thus we recommend users assess single-phenotype epistasis 
using functions in \texttt{cape} or in parallel analyses using tools such as \texttt{R/qtl} 
and \texttt{R/qtlbim} \citep{rqtlbim}.\\

\clearpage
 
 \section*{Tables of Functions}
 \addcontentsline{toc}{section}{Tables of Functions}
 \label{sec:Tables of Functions}
 
\begin{table}[h]\small
\caption*{\textbf{Basic \texttt{cape} Functions}}
\begin{tabular}{l p{10cm}}
\hline \hline
\multicolumn{2}{c}{\textbf{Reading and Writing Data}}\\
\hline
read.population & read in population data\\
writeVariantInfluences & write out a table of the variant influences to each other and the variant influences on traits.\\
\hline \hline
\multicolumn{2}{c}{\textbf{Data Manipulation}}\\
\hline
create.covar & assign a phenotypic value as a covariate in the genotype matrix\\
delete.pheno & delete the specified phenotypes from the phenotype matrix\\
get.covar & automatically select covariates for the pairwise scan based on a user-defined threshold of $\beta/\sigma$\\
norm.pheno & use rank Z normalization to normalize the phenotypes\\
select.eigentraits & select a subset of eigentraits to use in the analysis\\
select.pheno & select a subset of phenotypes to use in the analysis\\
select.by.chr & select a subset of chromosomes to use in the analysis\\
select.by.ind & subset the individuals in the population by either phenotypic or genotypic values\\
set.covar & manually set markers to be used as covariates for individual traits\\
set.pairscan.thresh & manually set the threshold for inclusion of markers in the pairscan\\
\hline \hline
\multicolumn{2}{c}{\textbf{Plotting Functions}}\\
\hline
plotNetwork & plot a network view of the significant influences of variants on each other and variants on traits\\
plotSinglescan & plot the results of the single-locus scan or markers selected for pairscan\\
plotPairscan & plot the results of the pair scan\\
plotSVD & plot the results of the singular value decomposition\\
plotVariantInfluences & plot the adjacency matrix showing the significant influences of variants on each other and variants on traits\\
\end{tabular}
\end{table}



\begin{sidewaystable}[h]\tiny
\caption*{\small{\textbf{Functions for cape Analysis}}}
\begin{tabular}{l P{3cm} P{3.5cm} P{1.5cm} P{3.2cm} P{2cm} P{4cm}}
\scriptsize{\textbf{Function}} & \scriptsize{\textbf{Behavior}} & \scriptsize{\textbf{Result Label}} & \scriptsize{\textbf{Structure}} & \scriptsize{\textbf{Internal elements}} & \scriptsize{\textbf{Each Row}} & \scriptsize{\textbf{Columns}}\\
\hline \hline
\nohyphens
calc.p & calculates empirical and corrected p~values of interactions using permutations & empirical.p & list of 2 matrices  &  m12, m21 & marker pair &  marker1, marker2, source, target, influence coefficient, influence $\sigma$, empirical p~value, adjusted p~value\\
 \hline
direct.influence & calculates the direct influence of each marker  pair on each trait and the p value of each influence & var.to.pheno.influence & list of  $N_t$$^*$ matrices & each matrix is named with its respective trait & marker pair & marker1, marker2, influence coefficients and standard errors for each marker on each trait\\
 &  & var.to.pheno.influence.perm$^{\dagger}$ & list of  $N_t$$^*$ matrices & each matrix is named with its respective trait & permutation of each marker pair & marker1, marker2, permuted influence coefficients and standard errors for each marker on each trait\\
 &  & var.to.pheno.test.stat$^{\dagger}$ & list of  $N_t$$^*$ matrices & each matrix is named with its respective trait & one marker in one marker pair context & marker, influence coefficient, $\sigma$, t~statistic, $\vert$t~statistic$\vert$\\
 &  & var.to.pheno.test.stat.perm$^{\dagger}$ & list of  $N_t$$^*$ matrices & each matrix is named with its respective trait & one permutation for each marker in one marker pair context & marker, influence coefficient, $\sigma$, t~statistic, $\vert$t~statistic$\vert$\\
 &  & max.var.to.pheno.influence & list of  $N_t$$^*$ matrices & each matrix is named with its respective trait & marker & marker, influence coefficient, $\sigma$, t~statistic, $\vert$t~statistic$\vert$, empirical p~value, adjusted p~value\\
 \hline
error.prop & propagates errors of coefficients calculated in the pair scan & var.to.var.influences or var.to.var.influences.perm & matrix &  & marker pair & marker1, marker2, m12, m21, m12~$\sigma$, m21~$\sigma$\\
\hline
get.eigentraits & performs singular value decomposition on phenotype matrix & ET & matrix &  & individual & one column for each eigentrait\\
\hline
maximum.influence & calculates the maximum influence that each marker has on each trait across all pair contexts & max.var.to.pheno.influence & list of  $N_t$$^*$ matrices & each matrix is named with its respective trait & marker pair & marker, influence coefficient, $\sigma$, t~statistic, $\vert$t~statistic$\vert$, empirical p~value, adjusted p~value\\
\hline
singlescan & performs all single-locus regressions for each marker on each trait & singlescan.results & list of  $N_t$$^*$ matrices & each matrix is named with its respective trait & marker & marker coefficient, $\sigma$, t~statistic, p~value\\
\hline
pairscan(n.perm $\geq 0$) & performs pair-wise marker regression for each trait & pairscan.results & list of  $N_t$$^*$ lists, each with 3 elements & twoD.effects & marker pair & marker labels and coefficients for each marker tested and each marker used as a covariate, including the interaction term\\
 &  &  &  & twoD.se & marker pair & marker labels and $\sigma$ for each coefficient in twoD.effects\\
 &  &  &  & model.covariance & marker pair & the model covariance matrix for each regression\\
 \hline
pairscan(n.perm $\geq 1$) & performs pair-wise marker regression and permutations for each trait & pairscan.perm & list of  $N_t$$^*$ lists, each with 3 elements & twoD.effects.perm & one permutation of marker pair & marker labels and coefficients for each marker tested and each marker used as a covariate, including the interaction term for all permutations\\
 &  &  &  & twoD.se.perm & one permutation of marker pair & marker labels and $\sigma$ for each coefficient in twoD.effects.perm\\
 &  &  &  & model.covariance.perm & one permutation of marker pair & the model covariance matrix for each permuted regression\\
 \hline
select.markers.for.pairscan & selects markers for the pair scan based on which markers are linearly independent of the other markers in the data set & geno.for.pairscan & matrix & & individual & one column for each independent marker selected for the pair scan\\
 &  & covar.for.pairscan &  matrix & & marker & each of the $N_t$$^*$ columns indicates which markers are to be used as covariates in the pair scan of each trait\\
 \hline
\end{tabular}
$^*$$N_t$ = Number of traits (phenotypes or eigentraits) scanned\\
$^{\dagger}$These objects are not added to the data object, but saved in \texttt{permutation.data.RData} if specified.
\end{sidewaystable}

\clearpage
 
\bibliographystyle{apalike}
\bibliography{cape}

\end{document}
